Uncertainties on the EFT coupling limits for direct dark matter detection experiments stemming from uncertainties of target properties

Direct detection experiments are still one of the most promising ways to unravel the nature of dark matter. To fully understand how well these experiments constrain the dark matter interactions with the Standard Model particles, all the uncertainties affecting the calculations must be known. It is especially critical now because direct detection experiments recently moved from placing limits only on the two elementary spin independent and spin dependent operators to the complete set of possible operators coupling dark matter and nuclei in nonrelativistic theory. In our work, we estimate the effect of nuclear configuration-interaction uncertainties on the exclusion bounds for one of the existing xenon-based experiments for all fifteen operators. We find that for operator number 13 the $\pm 1\sigma$ uncertainty on the coupling between the dark matter and nucleon can reach more than 50% for dark matter masses between 10 and 1000 GeV. In addition, we discuss how quantum computers can help to reduce this uncertainty and how the uncertainties are affected for couplings obtained for the nonrelativistic reductions of the relativistic interactions.


I. INTRODUCTION
The nature of dark matter (DM) has remained one of the biggest mysteries in physics.Over the past few decades, experiments have attempted to measure DM interactions with visible matter, both directly and indirectly, to no avail.Direct detection experiments for particle DM (see, e.g., Ref. [1] for a recent review), among which the currently largest and leading ones are XENONnT [2,3] and PandaX-4T [4], are generally built to search an area of parameter space consistent with a particular class of DM models, such as the Weakly Interacting Massive Particle (WIMP), see, e.g., Ref. [5] for a recent review.These experiments frequently look at only the leading order spin-independent and spin-dependent interactions [6][7][8][9], which may be suppressed by an asyet-unknown mechanism that prevents any detection of DM.Recently, there has been heightened interest to focus on a more general approach in the Galilean Effective Field Theory (EFT) [10,11] framework for setting constraints on DM-nuclei coupling and cross sections [12,13] for elastic scattering, as well as in Chiral Effective Field Theory for inelastic scattering [14].We note that while here we follow effective theories developed at the level of nucleon fields [10,11], other effective theories start at the level of quark and gluon fields [14][15][16][17][18][19][20][21][22].
Besides detectors built for the main goal of DM detec-tion, other low-energy, rare event search experiments can also be used in DM direct searches.Neutrinoless doublebeta decay (0νββ) is currently an extremely rare interaction of interest, and detectors such as CUORE designed for 0νββ also look for DM scattering events off their target material [23].CUORE uses TeO 2 crystals as a target; tellurium, like xenon (which we focus our analysis on in this paper), has isotopes that are currently modeled through nuclear configuration-interaction approaches because of the large number of nucleons (∼ 120-140).Thus, our calculations of WIMP-nucleon coupling uncertainties will be applicable to a wide range of sensitive, low-energy experiments.
In this paper, we show that a major source of uncertainty in the calculation of Galilean EFT DM-nucleus cross sections and the upper limits on DM couplings is nuclear models.Nuclei are highly complex manyparticle quantum systems and thus practically impossible to model or simulate to high accuracy without simplifying approximations on modern hardware.A common approach for heavier nuclei is the nuclear configurationinteraction in a shell model basis, where a frozen core of nucleons is surrounded by valence nucleons which occupy energy and angular momentum states in close analogy with atomic electrons [24].(Ab initio calculations of nuclei have been developed over recent years [25], but only very recently have these techniques started to be applied to heavy, open-shell nuclides such as xenon isotopes [26].)The target nuclear response to external probes such as DM scattering is then encoded in reduced matrix elements and a density matrix which contains the transition probabilities between different states of valence nucleon configurations.Still, this is a many-body system and thus challenging to calculate numerically.
In our work, we quantify the uncertainty in the up-, , FIG. 1.A schematic describing the general process to calculate upper limits on DM-nucleon EFT coupling coefficients and their uncertainties propagated from uncertainties in configuration-interaction calculations for a particular isotope.In our work, the two configuration-interaction Hamiltonians are GCN and JJ55 (described in Sec.II B).Our Monte Carlo procedure to construct Gaussian distributions for each one-body density matrix element (ρ f i J (a, b), see Eq. ( 5)) is described in detail in Sec.III.The differential DM scattering rates (dN/dEr, with N the event rate and Er the recoil energy, see Eq. ( 6)) are calculated using the dmscatter code which implements the intermediate calculational steps such as nuclear response function (W x,x ′ X , see Eq. ( 3)), differential cross sections (dσ/dEr, see Eq. ( 8)), DM halo model (f (v), see Eq. ( 7)), and DM response function R x,x ′ k (see, eg., App.C of Ref. [27], or Eq.38 in Ref. [10]).Then the energy integrated rates dN/dt were calculated adding the realistic XENON1T energy window, efficiency, and exposure to obtain the upper limits on the DM EFT coupling coefficients and their uncertainties as decribed in Sec.V.
per limits of the Galilean EFT DM-nucleus couplings by considering the uncertainty in configuration-interaction calculations of the target response.Figure 1 illustrates our process and provides a map of this paper.First, we choose as initial inputs two different appropriate configuration-interaction Hamiltonians, from which one generates the respective ground states and subsequently the relevant target one-body density matrices.From those computed density matrices we define a simple ensemble of density matrices, use a Monte Carlo (MC) procedure to sample matrices randomly from that ensemble, and characterize the resulting spread in the dependence of scattering event rates on recoil energy.Finally, we calculate the upper limits on the Galilean EFT WIMPnucleus couplings and their uncertainties coming from differences between the two model Hamiltonians.As discussed in Section III, this is a simplified approach to uncertainty quantification.
We find that the uncertainties in the upper limits of the coupling coefficients vary significantly between the different operators.The uncertainties for upper limits on couplings to some of the operators can be substantial, up to an order of magnitude for the ±1σ region.These results suggest that nuclear uncertainties need to be taken into account by experiments when placing upper limits on the DM-nucleus coupling coefficients for some operators.In addition, improvements in calculations of nuclear structure could significantly improve the accuracy of results emerging from dark matter detection experiments.
The remainder of this paper is organized as follows.
Sec. II presents the theoretical background used to calculate differential event rates for the 15 different possible nonrelativistic interaction channels between DM and nuclei.In Sec.III we present our methods to model the uncertainties coming from configuration-interaction calculations of the DM scattering rates.Sec.IV shows the results of the calculations for the ideal differential event rates and their uncertainties for different interaction channels.Then, in Sec.V we use our results together with knowledge of the properties of the XENONIT experiment to determine upper limits on the WIMP-nucleon coupling together with the propagated shell model uncertainties.In Sec.VI we compare the uncertainties arising from limits to our knowledge of the nuclear models to other sources of uncertainty.In addition, we show how the uncertainties on the coupling limits vary for the operators coming from nonrelativistic reductions of the relativistic interactions.We also discuss ways to improve the estimates of nuclear uncertainties as well as methods for reducing these uncertainties, including the possible utility of quantum computation in this context.In Sec.VII we conclude.review the nonrelativistic theory for DM-nuclei interactions, listing the derived operators for DM-nuclei elastic scattering.Next, in Sec.II B we show how the nuclear response functions are built from the nuclear density matrix elements, while in Sec.II C we present expressions for the differential cross-sections as a function of energy when the relevant nuclear response functions are known precisely.

A. Nonrelativistic Operators
The EFT calculations developed in Refs.[10,11], which assume Galilean invariance, lead to fifteen operators that are at most second-order in the exchanged momentum operators and are built as products of the following quantities: where mχ+m N is the dark matter-nucleon reduced mass with m χ (m N ) being the WIMP (nucleon) mass, ⃗ S χ is the WIMP spin and ⃗ S N is the nucleon spin, ⃗ q is the exchanged momentum, ⃗ v is the relative incoming velocity, and ⃗ v ⊥ • ⃗ q = 0 by momentum conservation.
These nonrelativistic EFT operators are [10,11]: where N = n,p for neutrons and protons, respectively.The total interaction Hamiltonian including the responses from all the operators is given by where the EFT coupling coefficients c x i , with i = 1, .., 15 being the operator index and x the nucleon type, are a priori unknown.In Sec.V we show how their values can be constrained by experimental non-detection of WIMPs with the data from the existing DM direct detectors.

B. Target response functions
The probability of a DM-nucleus interaction, i.e., the cross section, can be factorized into the dark matter response functions R x,x ′ k and the target response functions , where the index k = 1, .., 8 denotes the allowed combinations of electro-weak-theory operators.
The DM response functions, which group the operators derived from Galilean EFT theory (Eq.( 1)) for the corresponding nuclear response functions, can be found in App.C of Ref. [27], or Eq.38 in Ref. [10].
The target response functions are built from eight onebody electroweak multipole operators where Ψ i is the initial target wave function and Ψ f is the final target wave function, and J T is the angular momentum rank of the operator.The matrix elements of these operators, which can be constructed from Bessel spherical harmonics and vector harmonics [28], are calculated by summing over single particle orbitals from a one-body density matrix ρ: where a and b are the indices of the single-particle orbitals, and ⟨a||X x J T ||b⟩ are the reduced matrix elements of the operator X x J T , as obtained from the Wigner-Eckart theorem [29].
The target one-body density matrices are where ĉ † a and cb are the fermion creation and destruction operators [29].The wave functions and the subsequent density matrices used in this work were computed [27] using in configuration-interaction in a shell model basis via a high-performance code [30,31].
Since we assume the WIMP mass to be at least on the order of 10-1000 GeV and the WIMPs to be nonrelativistic, relatively little energy (O(100 keV)) is transferred to the nucleus.As such, we expect that the nucleus stays in the ground state throughout the interaction, i.e., |Ψ i ⟩ = |Ψ f ⟩.Thus, we are only interested in ground state-to-ground state transitions.
In our work, we use eigenstates and the resulting density matrices computed in a fixed Hilbert space using two sets of nuclear Hamiltonian matrix elements.The Hilbert space has a fixed core of 100 Sn and a valence space comprised of the 0g 7/2 -2s 1/2 -1d 3/2 -1d 5/2 -0h 11/2 orbitals.The two Hamiltonians are labeled as GCN [32,33] and JJ55 [34].The GCN shell model interaction was fit to experimental nuclear energies of nuclides in this mass region, starting from a G-matrix [35] constructed using the charge-dependent Bonn nucleon-nucleon (NN) potential version C, also known as the CD-Bonn potential [36].The JJ55 interaction is also obtained starting from a CD-Bonn G-matrix but was fitted to different experimental data from this mass region.Thus, the differences in the calculated target states represent the challenges in accurately modeling the nuclear many-body system.

C. Dark matter-nucleus event rates
The differential event rate dN/dE r (N is the number of events and E r is the recoil energy) is the convolution of the differential cross section dσ/dE r for the WIMPnucleus interaction and a DM halo model, integrated over DM velocity v: In the above equation, ε(E r ) is the detector efficiency, N T is the number of target nuclei in the detector, n χ is the local DM number density, and f (⃗ v), the DM velocity distribution in the rest frame of the detector, is obtained by boosting from the DM velocity distribution as seen in the rest frame of the Sun (the galactic-frame distribution).For this galactic-frame distribution, we use the simple halo model given by [37,38] f where Θ(v) is the Heaviside step function, v 0 = 220 km/s, N esc is a normalization factor, and the DM velocity cutoff v esc = 550 km/s, roughly the galactic escape velocity.The differential cross section, dσ/dE r in Eq. ( 6), for a DM-nucleus interaction for a particular isotope can be calculated as follows [27] where m T is the mass of the target nucleus (i.e., m T = Am N where A is the mass number of the nucleus and m N = 0.938 GeV/c 2 is the nucleon mass) and q = √ 2m T E r is the momentum transfer.The scattering transition probability [10,11,27] is given by where j T is the spin of the target nucleus and the sums are factored out into two parts, as mentioned in Sec.II B.
The first part contains the particle physics, i.e., the WIMP response function R x,x ′ i , and the second part contains the nuclear physics, i.e., the target response functions W x,x ′ X given in Eq. (3).

III. METHODS: MODELING UNCERTAINTIES IN DARK MATTER -TARGET EVENT RATES
In this section, we introduce our strategy for calculating the impact of uncertainty from the nuclear shell models on the DM-nucleus event rates.
To carry out uncertainty quantification (UQ) of DM detection is a daunting problem.A full UQ effort might involve hundreds or thousands of Monte Carlo (MC) samples of the input Hamiltonian matrix elements [39][40][41], the first box of Figure 1.The large configurationinteraction dimensions of the target isotopes, up to tens of billions, require thousands of hours of supercomputer time for each Hamiltonian sampled, an impractical methodology.Instead we choose to sample a model of the one-body density matrices, a procedure which may not fully reflect the uncertainties of the underlying Hamiltonian.Future work might use emulators, an increasingly popular approach in UQ [42], to speed up the sampling of the Hamiltonian, but emulators have yet to be fully applied to configuration-interaction calculations.
To estimate the impact of the uncertainties coming from configuration-interaction calculations on the upper limits of the WIMP-nucleon coupling, we use a simple MC procedure.As inputs to our MC simulations, we use two previously calculated nuclear structure data files containing reduced one-body density matrices calculated in the configuration-interaction shell model [27,30,31] for two different configuration-interaction Hamiltonians, GCN [32,33], and JJ55 [34], which are described in more detail in Sec.II B. This is represented in the flowchart in Fig. 1

by the box "Shell-model code."
To provide an approximate model of an ensemble of calculations, we calculate the mean and standard deviation of each density matrix element ρ f i J T (a, b) between the two interactions, i.e., for all combinations of initial and final orbitals a and b.With these averages and standard deviations, we construct a Gaussian distribution for each of the density matrix elements.Then, we randomly draw a value for each element from the defined distributions (This is in lieu of a MC sampling over an ensemble of Hamiltonians, which would be more fundamental but, as stated above, is not currently computationally tractable.) There is one caveat: the ground state-to-ground state density matrix elements for the J T = 0 transition for both protons and neutrons are each constrained by a sum rule in order to conserve particle number: where j a is the total angular momentum of the a-th valence orbital, J i is the total angular momentum of the state Ψ in Eq. ( 5), in our case always the ground state; Z valence and N valence are the number of valence protons and neutrons, respectively, and the bracket notation is We normalize the randomly generated J T = 0 density matrix element values such that the sum rules in Eqs.(10a) and (10b) are satisfied.There are no analogous sum rules for J T > 0, so it is possible to have unphysical values for those density matrices.Since our ensemble is defined by physical values, however, we nonetheless assume our calculations are sufficient to provide a realistic estimate of uncertainties.
The procedure described above can be briefly summarized with these three steps: 1. Construct a Gaussian distribution for each onebody density matrix element ρ f i J T (a, b) using matrices generated with at least two different nuclear shell models, 2. Draw randomized nuclear density matrix elements from created Gaussian distributions, 3. Normalize the new elements to satisfy particle number conservation, i.e., Eqs.(10a) and (10b).
The "Monte Carlo" section of the flowchart in Fig. 1 depicts the MC process we have just described.We use the generated random density matrix elements as inputs for the event rate calculations performed with the dmscatter code [27].We emphasize that this is a preliminary foray into uncertainty quantification (UQ) of the nuclear input into DM direct detection.Our results, given below, provide motivation for future investigations.

IV. RESULTS FOR IDEAL EVENT RATE SPECTRA
This section presents results for the ideal DM-nucleus differential event rates for xenon isotopes and their uncertainties calculated using the methods described in the previous section.
Figure 2 shows the 1σ and 3σ uncertainties in event rate spectra for WIMP scattering off 132 Xe and 129 Xe in the proton channel for a selection of EFT operators, assuming a WIMP mass of m χ = 150 GeV.Each of the event rate spectra subplots was generated with N = 1000 input files containing randomized nuclear density matrix elements produced from the Monte Carlo method described in the previous section.Several of the EFT operators yield significantly different event rates for the two shell models used.These large spreads translate into high uncertainties in the event rates.
Figure 2 emphasizes how the event rate spectra can differ significantly in shape and magnitude for various DM EFT operators.For example, the lower left and lower middle panels show the event rate spectra for operators O 4 -Eq.(1d) and O 6 -Eq.(1f) in 129 Xe.The former operator is the simple spin-dependent coupling where low energy recoils are favored (as in the case of coherent scattering), whereas the latter operator results in a peak in the spectrum due to an additional power of (⃗ q/m N ) 2 present in the cross section.
As mentioned, for many of the EFT operators, (see, e.g., the aforementioned O 4 , O 6 in Fig. 2,) the two configuration-interaction Hamiltonians GCN and JJ55 predict vastly different event rates, particularly in the low recoil energy regime; we hypothesize such differences are due to cancellations in Eq. ( 4).As such, the exclusion limits on the WIMP-nucleon coupling for these operators or the potential discovery fits will be sensitive to the choice of a given shell model interaction.
The mass of the WIMP affects the 'length' of the event rate spectra, i.e., where the event rate spectra tapers off to zero.As the WIMP mass increases, the WIMP carries more momentum and is thus able to cause collisions at higher nuclear recoil energies.This has implications for the discovery potential for DM detectors, as a 'longer' event rate spectrum would allow for more potential events if the detector is capable of registering events at higher nuclear recoil energies.On the other hand, a low WIMP mass would correlate to a 'short' event rate spectrum; if a detector is unable to probe scatterings at low recoil energies, it is possible that the majority of the WIMP-nuclei scatterings will go unseen, as the detector will be searching a recoil energy range that is mostly unpopulated.

V. SENSITIVITY ANALYSIS FOR XENON1T
So far there has been no conclusive detection of particle DM by a direct DM detection program.Therefore, only upper limits on WIMP coupling strengths to Standard Model particles have been set.
In this section, we describe how the process detailed in Sec.III can be used to calculate the uncertainty coming from nuclear shell models on the upper limits of WIMP couplings to nucleons.For the remainder of this work, we use the XENON1T experiment [7] as an example, although the uncertainty quantification procedure we will describe henceforth will work for any detector/element, provided that more than one nuclear density matrix data set is available.
XENON1T [7,43]  Recoil energy E r (keV) FIG. 2. Spread in the differential event rates dN/dEr for WIMPs interacting through DM nonrelativistic EFT operators O3, O12, and O15 with 132 Xe (top three plots), and O3, O12, O15, O4, O6, and O13 with 129 Xe (bottom six plots), assuming an effective exposure of 1 ton yr, ideal efficiency ε(Er) = 1, WIMP mass mχ = 150 GeV, and DM EFT coupling to protons c p i m 2 v = 0.25, where mv = 246.2GeV is the weak interaction mass scale.The plots demonstrate differences in the magnitude and shape of the scattering rates of WIMPs with Xe nuclei.Each panel is arbitrarily normalized such that the maximum of the GCN event rate (orange line) is 1/MeV, and the corresponding normalization factor C is provided in the top right of each plot.The actual ideal event rate is C multiplied by the event rate shown.
projection chamber DM direct search experiment located at the Laboratori Nazionali del Gran Sasso in Italy.The separate detection of the prompt scintillation in the liquid phase of the detector (S1 signal) and delayed scintillation in the gaseous phase of the detector (S2 signal) enabled a powerful discrimination between the electronic and nuclear recoils [14,44].As such, to good approximation the only backgrounds for the DM-nuclei interactions are those producing nuclear recoils only, which are close to zero in the 1 ton yr effective exposure of XENON1T [7,14].
To get the realistic upper limits on the DM EFT coupling coefficients we use the 1 ton yr effective exposure and detection efficiency ε(E r ) given by the curve in Ref. [14] used for the spin-independent analysis.We restrict our analysis to the region of interest in the recoil energy of E r = [4.9,54.4] keV.We calculate the total number of events in 1 ton yr exposure of XENON1T detector integrating the differential event rate Eq. ( 6) over the recoil energy window of the detector with the specified efficiency.
With the number of events per effective exposure, one can calculate the 90% confidence level (90% C.L.) upper limits on the EFT couplings of the DM particle.The 90% C.L. represents the smallest coupling, given a set of detector characteristics, at which a statistically significant number of WIMP-nuclei scatterings would be detected.Carrying this calculation out for varying WIMP mass m χ traces a curve in the c x i -m χ plane.The region of phase space above this curve represents all (c x i , m χ ) com-binations that through experimental non-detection have been excluded as realistic scenarios, while the region below represents combinations that have yet to be probed.To calculate these exclusion curves (upper limits), we assume that XENON1T is a zero background experiment, based on the successful discrimination of background signals in the S1 vs. S2 phase space plane [7,14].To model the detector composition, we use the abundance-weighted average of differential event rates over the six most abundant xenon isotopes ( 132 Xe, 129 Xe, 131 Xe, 134 Xe, 136 Xe, and 130 Xe, in that order).The relative abundances of these isotopes by mass are respectively 26.9%, 26.4%, 21.2%, 10.4%, 8.86%, and 4.07%; when summed, these six isotopes constitute slightly more than 97% of the mass of the XENON1T detector.There is a different amount of uncertainty that each isotope contributes, which is discussed further in Sec.VI.
The event rate dN/dt is linearly dependent on the dark matter response functions R x,x ′ X through the cross section.Further, if we assume isospin is conserved x = x ′ (i.e., we only look at p→p and n→n interactions as is the case in the ground state-to-ground state interaction), the R x,x ′ X depend quadratically on the EFT coefficients.If only a single EFT coefficient c x i is nonzero then the total event rate for that coupling is a function of (c x i ) 2 .Thus, to find the minimal coupling needed to produce a zero background event rate dN min /dt required for a statistically significant WIMP-nucleus scattering detection, we only need to run the above event rate procedure once for a test coupling c x i,0 .Assuming the number of WIMP-nucleon scattering events follows the Poisson distribution, in order to achieve a 90% confidence of detection we need an event rate of at least dN min /dt = 2.3 events ton −1 yr −1 in the zero background case.The chosen test coupling c x i,0 yields an event rate of dN 0 /dt, which we then use to scale the EFT coefficient to find the minimal required coupling for detection, c x i,min : By calculating the minimal couplings for WIMP masses in the range 10 GeV ≤ m χ ≤ 1000 GeV, we create exclusion curves for each EFT coefficient as a function of m χ for repeated calculations with randomized nuclear density matrices.We ran N trial = 1000 MC trials for most of EFT coefficients for each m χ 1 .Thus, we can quantify the uncertainty on the coupling limits that derive from the uncertainties in the nuclear models from the spread in the N trial exclusion curves.
The distributions of the exclusion curve values at each m χ are in general not Gaussian despite the nuclear matrix elements in the MC simulation being drawn from 1 For O 13 , the distribution of exclusion curves was highly non-Gaussian, and as such we increased N trials to 10 4 in order to obtain a more accurate probability distribution.
normal distributions.To calculate the 1σ bands, we use the Feldman-Cousins procedure described in Ref. [45].We find the interval in coefficient values [a, b] such that the integral of the histogram and H(a; m χ ) = H(b; m χ ). Figure 7 in App.B shows an example of such a 1σ band for c p 13 at m χ = 12 GeV.Figure 3 shows the calculated examples of the upper limits on the DM-nuclei coupling coefficients for protons (upper panels) and neutrons (lower panels) as a function of WIMP mass and their uncertainty for selected DM EFT operators calculated for the XENON1T detector.The remaining combinations that exhibit sensitivity to the shell-model uncertainties are presented in App.A for completeness.
The WIMP-nuclei upper coupling limit is noticeably more sensitive to the density matrices when coupling through operator O 13 in Eq. (1m) than other EFT DM operators for both protons and neutrons, corresponding to potentially larger uncertainty.Notably, for most of the neutron EFT operators, the exclusion curves generated using our MC process (see Sec. III) are normally distributed.However, the distribution for O 13 is not normally distributed and is skewed toward lower couplings.Operators O 9 and O 14 for neutrons and O 4 , O 6 , O 9 , O 10 , and O 14 for protons also exhibit relatively greater sensitivity for both proton and neutron couplings.

VI. DISCUSSION
Our results demonstrate that the need to implement nuclear models that involve approximations can lead to substantial uncertainties in the interpretation of the results of dark matter detection experiments.
Figure 4 shows the upper limits on the nonrelativistic proton-WIMP coupling coefficient values and their relative uncertainties for 129 Xe, 131 Xe and 132 Xe assuming m χ = 70 GeV, and the total fiducial mass of XENON1T is comprised of each isotope, respectively.We show here the results only for the most abundant even-even isotope of xenon because all four even-even isotopes had nearly identical sensitivities and uncertainties for each EFT coefficient and did not probe operators O 4 , O 6 , O 7 , O 9 , O 10 , O 13 , or O 14 .Standard selection rules for parity and angular momentum cause the matrix elements of these operators to vanish for J i = 0 ground states.
The two even-odd isotopes, 129 Xe and 131 Xe, have J i > 0 ground states and thus probed all operators while producing more uncertain upper limits in most cases for the operators present for even-even isotopes.
Although the dark matter-nucleon couplings were derived from Galilean-invariant effective field theory [11], one can also start from a relativistic, Lorentz-invariant WIMP mass m χ (GeV/c 2 ) 90% CL upper limit FIG. 3. The 90% C.L. upper limits on the nonrelativistic EFT coupling coefficients (as detailed in Sec.V) for select nucleon and operator combinations along with 1σ and 3σ uncertainty bands from nuclear uncertainties calculated using XENON1T experiment as a model for a general DM direct detection experiment.The top three plots are exclusion plots for the EFT couplings to protons, while the bottom three plots are exclusion plots for the EFT couplings to neutrons.We use a common vertical scale for all of the plots to show which operators are more constrained than others.Note that the vertical scale obscures the non-zero sensitivity of O4 Eq. (1d) to the variations in nuclear density matrices (bottom left plot).
Lagrangian, including those derived from chiral effective field theory, and carry out a nonrelativistic reduction [11,16,[19][20][21][46][47][48] of the WIMP-nucleon interaction terms.This leads to specific linear combina-tions of the nonrelativistic operators, with possible additional dependencies upon the momentum transfer q and WIMP mass m χ .We considered the seven cases given in Ref. [48]: where χ is the WIMP four-spinor (and thus one must assume the WIMP is spin 1/2) and N the nucleon spinor; m M defines the scale of the effective theory, often chosen to be either the weak interaction scale, m v ≈ 246 GeV [11], or the nucleon mass [48].
The linear combination of the Galilean EFT operators allows for potential interference between operators when calculating event rates.Coupled with the momentum transfer and m χ dependence, uncertainties in minimal relativistic couplings C x i for the relativistic interaction terms could have a non-trivial relation to the minimal coupling uncertainties of the Galilean EFT operators.Thus, we repeated the analysis of the uncertainty in minimal coupling coefficients (see Secs.III -V and Fig. 4) for the seven relativistic operators L int i above.To do this, we modified the dmscatter code [27] to specify these relativistic couplings.
Figure 5 shows the 90% confidence level minimal coupling coefficients and their relative uncertainties for the relativistic interaction couplings.As expected, terms that are combinations of Galilean EFT operators with negligible uncertainties (i.e., L int 5 , L int 9 , L int 13 , and L int 17 ) have negligible uncertainties derived from the nuclear structure uncertainties.In general, the uncertainties do not reflect significant interference between the nonrelativistic operators; the one exception is interaction term 10, which has a relative uncertainty in 131 Xe that exceeds the uncertainty, by a factor of roughly 3, in either of the Galilean EFT operator minimal couplings that it depends on, namely O 4 and O 6 .Therefore, our conclu-sions are not very sensitive to the choice of relativistic or nonrelativistic formalism.
Figures 4 and 5 suggest that terms containing nuclear spin ( ⃗ S N in Eqns.(1) and nucleon axial vector or tensor terms in Eqns.( 13)) have the largest uncertainties stemming from nuclear model uncertainties.In Appendix C we compare energies and electromagnetic observables from our two models to experimental results for pertinent xenon isotopes.Among the observables are magnetic dipole moments, which include an effective coupling to nucleon spins.We get reasonable agreement, with neither model systematically better than the other; is most cases the calculated electromagnetic observables are within experimental error.This bolsters confidence in our ability to model the nuclear response to WIMP scattering, including couplings dominated by nucleon spin, as well as to estimate theoretical uncertainties.
Uncertainties from nuclear modeling are obviously not the only source of error in dark matter detection experiments.For example, recoil rates can be heavily sup- pressed by isospin-violating DM interactions in DM direct detection experiments [49,50].Similarly, Ref. [51] showed that non-parametric uncertainties in the dark matter halo velocity distribution could weaken exclusion limits in XENON1T by as much as two orders of magnitude, although this assumes m χ ≤ 60 GeV and anisotropy in the dark matter halo.Above m χ = 60 GeV, uncertainties have a mild effect on the limits.Therefore, there is a reasonable parameter space in which the nuclear uncertainties studied in this paper could dominate the overall uncertainty of the experiments.

A. Strategies for improving the calculations
The methods used to obtain the results presented here could be improved both by the use of more sophisticated but computationally intensive methods to estimate the uncertainties and by developing methods that enable the uncertainties themselves to be reduced.
Regarding our estimates of the uncertainties, it is worth emphasizing that we have modeled the uncertainty in the nuclear shell-model by creating variations in the reduced one-body density matrix elements, as shown in the center column of Fig. 1.A more precise uncertainty analysis would begin by estimating the uncertainties in the nuclear Hamiltonian matrix elements [40] (top-left box of Fig. 1), and propagating this uncertainty [41] through the same flowchart from there, ultimately to the minimal EFT couplings again.The downside to this alternate method is that the configuration-interaction calculations used to reduce the many-body system to the more tractable one-body density matrices (see left column of Fig. 1) are computationally intensive as well as time-consuming.One potential alternative is the use of emulators [42], but the application to such shell-model calculations has yet to be performed.
Much of the challenge of modeling the nuclear response to DM is the daunting dimensions, to 10 10 or beyond, of even a truncated Hilbert space.Improving the uncer-tainty estimates and reducing the uncertainties of the calculations themselves are both areas in which quantum computing has the potential to excel.We note, however, the large number of nucleons in the relevant nuclei and the computational depth of the necessary computations are such that achieving higher accuracies than classical methods would require substantial improvements in quantum computing hardware and error correction.In any case, the preliminary analysis that we have presented here can highlight which isotope and EFT operator combinations exhibit nontrivial uncertainty and therefore are attractive targets for further investigation.

VII. CONCLUSIONS
Understanding the nature of dark matter is critical to progress in physics, as it could shed light on some of the most fundamental questions in cosmology and particle physics.While there has been significant progress in direct dark matter detection experiments, so far only gravitational interactions of dark matter have been observed.Direct detection experiments aim to observe the interactions of particle dark matter with ordinary matter, which could provide valuable information about its properties, such as mass, type of the interaction, and interaction strength.This information can help in determining the underlying theory of dark matter, which could have significant implications for our understanding of the universe at large.
We have presented a preliminary uncertainty analysis of the WIMP-nucleon coupling exclusion limits stemming from uncertainties in nuclear shell model calculations.We have applied this 'first-order' analysis to study the sensitivity of all operators in the Galilean effective field theory, and identified the EFT operators (couplings) which are particularly affected by current limitations of nuclear shell model calculations.The analysis presented here could also be easily applied to direct detection experiments using different elements in the detector, such as argon (DarkSide-50 [12]) or tellurium (CUORE [23]).We also calculated the change in the uncertainties on the coupling limits for the operators obtained by performing nonrelativistic reductions of the relativistic interactions.
Our analysis reveals that the degree of uncertainty in the upper limits of the coupling coefficients varies considerably across the different Galilean dark matter Effective Field Theory operators.Specifically, we observe that the uncertainties for the upper limits on couplings to certain operators can be substantial, with the ±1σ range spanning over an order of magnitude in the limit of low dark matter masses below approximately 20 GeV for the operator O 13 .These findings underscore the importance of accounting for nuclear uncertainties when determining the upper limits on dark matter-nucleus coupling coefficients in direct dark matter experiments.Moreover, the advancements in the precision of nuclear structure calculations have the potential to greatly enhance the ac-curacy of results obtained from dark matter detection experiments.FIG. 6.The C.L. upper limits on the nonrelativistic EFT coupling coefficients (as detailed in Sec.V) for select proton and operator combinations along with 1σ and 3σ uncertainty bands from nuclear uncertainties calculated using XENON1T experiment as a model for a general DM direct detection experiment.Note that the vertical scale obscures the non-zero sensitivity of O3 Eq. (1c), O12 Eq. (1l) and O15 Eq. (1c) to the variations in nuclear density matrices.

Appendix B: Probability distributions for couplings
In this appendix we include an example plot illustrating the Feldman-Cousins [45] procedure for finding the uncertainties especially useful for non-symmetric distributions.
Figure 7 shows the distribution of the calculated 90% C.L. upper limits on the coupling coefficient of operator O 13 for protons calculated for N trials = 10 4 MC generated one-body density matrix files with the procedure described in Sec.III.The obtained distribution of the coupling coefficients is highly non-normally distributed although each of the one-body density matrix elements distributions were constructed with Gaussians.

Appendix C: Tests of nuclear models
Our Monte Carlo sampling is built upon two input models, described in Section II B. Both models have been previously compared against experimental data.The original JJ55 paper [34] found good agreement with excitation spectra of selected Sn, Sb, and Te isotopes as well as magnetic dipole (M1) moments of selected states of Sn, Te, Xe, and Ba isotopes.The original GCN5082 papers [33] found good agreement with experimental spectra of light Xe isotopes.Several subsequent papers compared properties using both models: spectrum, transition strengths, and g factors for 130 Te [60]; and the energy spectra (levels) of 130 Ba [61] and of 131 Xe [62].Using JJ55, comparisons to experiment have been made of electric quadrupole (E2) and magnetic dipole (M1) transitions in Te isotopes [63,64] and in 131 I and 132 Xe [65].Finally, we note that GCN5082 has been used previously for calculations of dark matter cross sections [66].
Here we provide additional comparison of observables of those two models to experimental results.Specifically we compare against the measured low-lying excitation spectra, electric quadrupole gamma decays, and static magnetic dipole moments.Following common practice for valence shell calculations [67], we use effective charges and g-factors, with values taken from the literature, and do not use any explicit current corrections.[34] and GCN5082 [33].All energies are in MeV.Here the subscript n denotes the order in the spectrum for a given J π : thus 2 + 1 is the lowest 2 + level and 2 + 2 is the next lowest.All Xe isotopes with even A have an empirical J π = 0 + ground state, while 129 Xe has a 1/2 + ground state and 131 Xe has a 3/2 + ground state.Note that the JJ55 model puts the 1/2 + 1 level below the 3/2 + 1 level, at odds with experiment, but the discrepancy is within the typical theory uncertainty, about 126-150 keV, of similar configuration-interaction calculations [40,59].The experimental errors on the energies are sub-keV.
In Table II we compare experimental reduced electric quadrupole matrix elements, or B(E2)s, against calculated values.For our calculations, we used the standard prescriptions [24,67], working in a single-particle harmonic oscillator basis with an oscillator length parameter of 2.295 fm, and adopted the effective charges of 1.86 e and 0.65 e for valence protons and neutrons, respectively [63].We only considered even mass numbers A to avoid the complication of mixing E2 and M1 transitions.
Finally, in Table III, we compare experimental static magnetic dipole moments, or µ, against calculated values.Following [34], we quench the free spin g factors by 0.7.
Broadly speaking, both models give similar agreement with experiment, justifying using them as an equally weighted basis for our calculations.
was a 3.2 ton dual-phase xenon time

FIG. 4 .
FIG.4.The 90% C.L. upper limits on the proton-WIMP nonrelativistic EFT coupling coefficients (top) and their relative uncertainties (bottom) for 129 Xe, 131 Xe, and 132 Xe, assuming a WIMP mass of mχ = 70 GeV.Not shown are the three other even-even isotopes as they have nearly identical sensitivities to 132 Xe for all EFT channels.Note that xenon isotopes with even mass number are not sensitive to operators O4, O6, O7, O9, O10, O13, and O14.

FIG. 5 .
FIG.5.The 90% C.L. upper limits on the proton-WIMP relativistic couplings (top) and their relative uncertainties (bottom) for 129 Xe, 131 Xe, and 132 Xe, assuming a WIMP mass of mχ = 69.5 GeV.Definitions for the relativistic couplings in terms of the nonrelativistic EFT operators can be found in Eq.(13).Note that even-even xenon isotopes are not sensitive to operators L int 7 , L int 10 , and L int 15 .

Distribution of c p 13 m 2 vFIG. 7 .
FIG. 7.Example distribution of the 90% C.L. upper limit on c p 13 for N trials = 10 4 MC density matrices at a WIMP mass of mχ = 12 GeV.The shaded green region represents the 1σ band and contains 68% of the distribution as found from the Feldman-Cousins procedure as explained in Sec.V.The 3σ band, containing 99.7% of the distribution, is approximately represented by the support of the histogram.The obtained distribution is highly non-symmetric.

TABLE I .
Table I compares experimental excitation energies of low-lying states against the values calculated in our two models.Note that the typical theory error in shell model Excitation energies in Xe isotopes, comparing experiment against two models, JJ55