Nuclear structure factors for general spin-independent WIMP-nucleus scattering

We present nuclear structure factors that describe the generalized spin-independent coupling of weakly interacting massive particles (WIMPs) to nuclei. Our results are based on state-of-the-art nuclear structure calculations using the large-scale nuclear shell model. Starting from quark- and gluon-level operators, we consider all possible coherently enhanced couplings of spin-1/2 and spin-0 WIMPs to one and two nucleons up to third order in chiral effective field theory. This includes a comprehensive discussion of the structure factors corresponding to the leading two-nucleon currents covering, for the first time, the contribution of spin-2 operators. We provide results for the most relevant nuclear targets considered in present and planned dark matter direct detection experiments: fluorine, silicon, argon, and germanium, complementing our previous work on xenon. All results are also publicly available in a Python notebook.


I. INTRODUCTION
Astrophysical observations have established that more than three quarters of the matter content of the universe are composed of dark matter. The nature of dark matter, however, remains elusive, and its very existence is one of the most compelling pieces of evidence for physics beyond the Standard Model of particle physics. A key step to unveil the composition of dark matter would be its direct detection in the laboratory [1]. Such endeavor is led by international collaborations that use atomic nuclei as targets, in their aim to detect the nuclear recoils resulting from the scattering of dark matter particles [2][3][4][5][6][7][8][9][10][11][12][13]. In spite of the very high sensitivities achieved by reduced backgrounds combined with extended exposures, there has been no conclusive evidence to date for the direct detection of dark matter. Next generations of experiments plan to push the frontier for dark matter direct detection by several orders of magnitude [14][15][16][17], until the background from coherent neutrino-nucleus scattering [18] becomes dominant.
Direct detection experiments are motivated by extensions of the Standard Model that propose dark matter candidates interacting with quarks and gluons, the Standard Model fields that ultimately form atomic nuclei. Prominent candidates are weakly interacting massive particles (WIMPs) [19]. Because the WIMPs forming dark matter would be non-relativistic (NR), their scattering off atomic nuclei would transfer energies and momenta much smaller than the nuclear or nucleon masses. As a consequence nucleons and nuclei, instead of quarks * E-mail: mhofer@uw.edu † E-mail: pklos@theorie.ikp.physik.tu-darmstadt.de ‡ E-mail: menendez@cns.s.u-tokyo.ac.jp § E-mail: schwenk@physik.tu-darmstadt.de and gluons, become the relevant degrees of freedom. The interpretation of the present experimental limits and future values of WIMP-nucleus cross sections, therefore, naturally depends on the nuclear physics aspects of the scattering. This information is encoded in the so-called nuclear structure factors [20]. While for some simple cases a phenomenological prescription can be a good approximation to the structure factor [21,22], in general a good description of the nucleus, obtained with a dedicated nuclear many-body calculation, is needed.
In the absence of experimental data on the spin [23,24], momentum-transfer [25,26], or isospin dependence, the character of the WIMP-nucleus interaction is unknown. Nevertheless standard analyses usually assume so-called spin-independent (SI) interactions, which receive the coherent contribution of all nucleons in the nucleus. However, additional interactions are possible, and could be dominant if the SI coherence is compensated by suppressed values of the corresponding WIMP-nucleon couplings. Some examples analyzed experimentally are the so-called spin-dependent (SD) [3,[27][28][29][30][31] or momentumtransfer-dependent interactions [32].
In order to organize different possible WIMP-nucleus interactions, two alternative schemes have been proposed recently. On the one hand, a non-relativistic effective field theory (NREFT) [33][34][35] based on the lowest-order operators that can describe the coupling of a WIMP to a nucleon. The NREFT approach proposes a set of one-body operators, considered in recent experimental analyses [36][37][38][39]. On the other hand, an organization based on chiral effective field theory (ChEFT) [40][41][42], a low-energy effective theory of quantum chromodynamics (QCD) that preserves the QCD symmetries, in particular capturing the important role played by pions at low energies. The ChEFT approach can be mapped onto the single-nucleon couplings of NREFT [43], implying certain interdependencies for the latter. In addition, ChEFT predicts consistent couplings of WIMPs to two nucle-ons [44][45][46][47][48][49][50][51], reflecting that nucleons are strongly interacting in nuclei. Such contributions occur, e.g., when the WIMP couples to a virtual pion exchanged between the two nucleons, an effect recently constrained for the first time by direct-detection experiments [52], for the case of a scalar WIMP-quark interaction. Similar couplings to two nucleons are important in electromagnetic and weak transitions of atomic nuclei [53][54][55][56]. Furthermore, ChEFT provides a power counting that suggests a hierarchy, guided by QCD, for the expected importance of the different NREFT operators [43]. This hierarchy is only tentative, because the couplings describing the interaction of the WIMP to the Standard Model fields are not known. A related ChEFT approach limited to WIMP interactions with one nucleon has been proposed in Refs. [57,58].
In this work we follow Ref. [49] and combine the ChEFT framework with large-scale shell model nuclear many-body calculations to calculate the leading nuclear structure factors that exhibit the coherent contribution of several nucleons in the nucleus. We call these generalized SI interactions. We consider both WIMP couplings to one and two nucleons, with special emphasis on the two-nucleon contributions related to the diagonal and non-diagonal parts of the energy-momentum tensor. We note that at this point our structure factors are not fully consistent, in the sense that our many-body calculations are based on phenomenological interactions instead of ChEFT. Such consistent studies are presently only available for few-nucleon systems [59]. In addition, nuclear response functions based on nuclear states obtained using ChEFT interactions are available for light nuclei [60]. However, such very light isotopes are not used in leading direct detection experiments. While Ref. [49] was limited to a xenon target, here we study all the stable isotopes of fluorine, silicon, argon, and germanium, the nuclear targets used and considered in present and future direct detection experiments. In addition, we provide a Python notebook to facilitate the use of our structure factors for both theorists and experimentalists.
The rest of the article is organized as follows. In Sec. II we describe the formalism and discuss which terms we include in the cross section, focusing on a spin-1/2 WIMP. Section III introduces the nuclear structure many-body calculations that describe the nuclear targets used in direct-detection experiments. The results for our calculated structure factors are given in Sec. IV. We distinguish between WIMP couplings to one nucleon, discussed in Sec. IV A, and the coupling to two nucleons, the subject of Sec. IV B. The most important new aspects of our implementation of two-body effects are highlighted in Sec. V. We conclude with a summary of the main findings of this work in Sec. VI. Details of the matching to NREFT, the nucleon matrix elements, matching for a spin-0 WIMP, and the calculation of two-body interactions are provided in Appendices A-C. All our results are also available in the form of a Python notebook, with a brief User's guide in Appendix D.

II. FORMALISM
Based on Ref. [49] we consider the following cross section for the generalized SI WIMP-nucleus scattering: where q = |q| is the three-momentum transfer between the WIMP and the nucleus, and m χ , m N , and m N denote the WIMP, nucleon, and nuclear masses, respectively. The relative velocity of the WIMP is v = |v|, [34] is the velocity of the WIMP with respect to the nucleus, with reduced mass µ N = m N m χ /(m N + m χ ). Each term in Eq. (1) contains a structure factor F and a coupling c. The structure factors F encode the nuclear structure aspects of the scattering, and are the main subject of this work. The couplings c are a convolution of the Wilson coefficients, which describe the fundamental interaction of the WIMPs with the quarks and gluons, and the hadronic matrix elements. They therefore depend on the particular scenario for physics beyond the Standard Model considered. Finally, the kinematic factors with the reduced mass µ N = m N m χ /(m N + m χ ), appear in the contribution from the subleading NREFT operators O 5,8,11 , defined in Appendix A.

A. Effective Lagrangian and couplings
To be definite, we consider the case of a spin-1/2, Standard-Model singlet χ throughout the main text, but to demonstrate that the decomposition in Eq. (1) applies in full generality we also provide the matching relations for a spin-0 WIMP, see Appendix B. For spin-1/2, the relevant terms in the effective Lagrangian [61][62][63] are where we have only listed operators that can lead to coherently enhanced responses or feature in the standard SD interaction, and χ we already integrated out the heavy quarks [64], whose effect is absorbed in while elsewhere the sum runs, in principle, over all quark flavors q. C S g andC S g are the original coefficients of the gluon operatorχχ α s G a µν G µν a andχiγ 5 χ α s G a µν G µν a , respectively, rewritten in terms of the trace of the energymomentum tensor: Moreover, we introduced its traceless componentsθ µν = θ µν q +θ µν g in the context of the dimension-8 spin-2 contribution:θ In detail, the coefficients in Eq. (1) are given by (nucleons N encompass protons and neutrons, N = p or n) where ζ = 1(2) for a Dirac (Majorana) particle. 1 The various couplings refer to combinations of Wilson coefficients from Eq. (3) and nucleon matrix elements The explicit expressions, together with the matching onto NREFT, are collected in Appendix A.

B. Pion matrix elements
Here, we highlight the couplings to the pion in Eq. (7): In this case, the relevant matrix elements are the following: first, the scalar couplings to u-and d-quarks are given by where [65] (and f π s ≈ 0). Second, the coupling to the trace anomaly is only sensitive to the sum f π u + f π d = 1, so that no new coupling appears in f θ π . The spin-2 couplings are related to moments of pion parton distribution functions (PDFs) q π (x): and similarly for the gluon, subject to the constraint q f (2) q,π + f (2) g,π = 1.

C. Dipole operators
For a Dirac WIMP, a possible extension beyond Eq. (3) concerns dark matter candidates with a non-vanishing dipole moment, corresponding to the effective dimension-5 Lagrangian [71,72] where due to σ µν γ 5 = i 2 µναβ σ αβ the second term involving the dual field strength tensorF µν = 1 2 µναβ F αβ is equivalent to −χσ µν iγ 5 χF µν . These operators produce long-range tree-level interactions via the exchange 2 Note that all spin-2 couplings are scale-dependent quantities, so that, at a scale of µ = 2 GeV in MS, the value from Ref. [68] actually becomes x π v = 0.256(13) [69]. This has been taken into account in Eq. (14). of a photon, which leads to the NR one-body amplitudes Here, F N 1/2 (t) are the Dirac/Pauli form factors of the nucleon with the full dependence on the relativistic momentum transfer t = −q 2 (up to relativistic corrections). The corresponding extension of Eq. (1) is straightforward, with the photon poles introducing new terms, besides additional contributions to some of the existing ones. Formally, this amounts to (q-dependent) terms in the matching relations in Eq. (7): where we have already taken ζ = 1 due to the absence of tensor currents for Majorana particles. The Dirac-formfactor radii are related to the Sachs ones by In the spirit of the present study, Eq. (18) neglects the non-coherent contributions from O 4 and O 6 . Further, there are many more dimension-7 operators involving the (electroweak) field strength tensors [73,74], and matching relations similar to Eqs. (17) and (18) could be extended accordingly.
In the form (17) the amplitudes automatically include radius corrections, subsumed in the full electromagnetic form factors. In addition, Eq. (16) could in principle produce new two-body currents. In the end, such terms take a similar form as the axial-vector-vector two-body currents identified in Ref. [43], and are only suppressed by two chiral orders compared to the leading pieces of Eq. (17). However, as argued in Ref. [49], after summation over spins their isospin structure [τ 1 × τ 2 ] 3 leaves only an isovector coherent enhancement suppressed by (N − Z)/A with respect to the scalar two-body current. We will continue to neglect such effects in the remainder of this work.

D. Scaling of operators
The global picture that arises in this way is summarized in Table I. For each nuclear structure factor we list: Structure factor QCD operators Chiral scaling NR operators Overall scaling Interference with O1 Likewise, the table only shows coherently enhanced contributions for operators with velocity v ⊥ , which in nuclei produce both a term that behaves as q/mN and a remainder determined by the target velocity v ⊥ T ∼ 10 −3 (see Appendix A for more details). For comparison the last rows for Sij show operators that lead to SD interactions. In the chiral counting we have Mπ = O(p), v ⊥ = O(p 2 ), with relativistic corrections counted as ∂/mN = O(p 2 ) and ∂/mχ = O(p 2 ). Two-body structure factors, Fπ and F b , cannot be matched onto single-nucleon NREFT operators, and the quasi-coherence of F Φ is characterized by ξ ∼ 0.2. Finally, the last column indicates whether each entry interferes with the leading O1 operator. the relativistic operators that contribute, the corresponding NREFT operator if applicable, the chiral scaling as well as the overall scaling including coherence, and finally whether or not the respective contribution interferes with the leading O 1 operator. Table I reflects the strategy underlying the present work, in that it includes all contributions, either twobody currents or subleading one-body operators, that appear up to third order in the chiral expansion and receive some form of coherent enhancement. Selected coherent chiral fourth-order contributions are shown as well, mainly because in many cases this is the order when some of the coherent NREFT operators first enter. For comparison, the table also includes the leading relativistic operators that produce the NR expansion related to the standard SD interactions, even if these are not coherent.

III. NUCLEAR STRUCTURE CALCULATIONS
The calculation of the nuclear structure factors requires a many-body approach that describes the ground states of the target nuclei considered. As in previous works [22,23,46,47,49] we use the nuclear shell model, one of the most successful many-body approaches in medium-mass and heavy nuclei [75]. For all calculations we have used the shell model code ANTOINE [75,76].
The shell model is based on the solution of the quantum many-body problem in a reduced configuration space where the Schrödinger equation for the nuclear ground and low-energy excited states can be solved exactly. We highlight two important aspects. First, the configuration space used in the calculation needs to capture the nuclear structure properties relevant for the process of interest. The limitation of using a restricted configuration space stems from the difficulty to solve the nuclear many-body problem in a nontruncated space for heavier nuclei. For nuclear targets used in direct detection experiments, cal-culations of structure factors without such truncations exist up to 4 He [51,59,60] and could be performed in the near future up to 40 Ca. Second, calculations must use an effective interaction appropriate for such a configuration space.
Until very recently, nuclear shell-model calculations relied on phenomenological effective interactions. In spite of being derived from nucleon-nucleon (N N ) scattering data, these interactions have to be adjusted phenomenologically, mostly the part that describes the singleparticle aspects of the nuclear interaction, referred to as monopole part, in order to achieve a better agreement with the nuclear structure of heavier nuclei. Progress in nuclear theory has improved this picture by including in the starting point of the derivation of effective interactions, in addition to N N interactions, also nuclear forces between three nucleons, 3N forces. Three-nucleon forces are the analog of two-body currents in the coupling to external probes, such as WIMPs, to nucleons. Consistent N N and 3N interactions can be derived in the ChEFT framework [40][41][42]. Starting from these ChEFT N N and 3N forces it is possible to derive effective interactions that, without further phenomenological adjustments, reproduce nuclear spectroscopy rather well in nuclei with nucleon number A comparable to direct-detection nuclear targets [77][78][79]. First studies have just started to extend these techniques to study electromagnetic and weak transitions in medium-mass nuclei [80][81][82]. In addition, the effective-theory character of ChEFT, combined with the consistency between nuclear forces and currents, provides a framework to quantify nuclear-structure uncertainties [83][84][85]. In this work, however, we follow the standard approach and use phenomenological effective interactions. This prevents the assessment of reliable nuclear-structure uncertainties, which will become possible in the future.
The configuration spaces and effective interactions we have used are described below. All single-particle orbitals belong to a three-dimensional harmonic oscillator basis nl j , where n is the principal quantum number, and l, j denote the orbital and total angular momentum.
The lightest target we have studied has only one stable isotope, 19 F, which we already considered in Ref. [47]. Here we use the same USDB effective interaction and the 0d 5/2 , 1s 1/2 , and 0d 3/2 single-particle orbitals, for both neutrons and protons. This configuration space is known as the sd shell. Reference [47] used the same configuration space and effective interaction for 29 Si, and here we extend the study to the other two stable silicon isotopes, 28,30 Si.
For argon, we study only 40 Ar. The configuration space we consider is significantly larger, and comprises seven single-particle orbitals for neutrons and protons, the three sd-shell orbitals plus the 0f 7/2 , 1p 3/2 , 0f 5/2 , and 1p 3/2 , where the latter comprise the so-called pf shell. We use the SDPF.SM effective interaction, which describes well the electromagnetic properties of ground states and the coexistence of spherical and deformed  states in this mass region [86,87]. In order to make the diagonalizations in the configuration space feasible, we need to truncate our many-body calculations, by keeping the 0d 5/2 orbital filled with nucleons, and restricting the number of excitations from sd-shell to pf -shell orbitals to 8. These are similar truncations to those in Ref. [86], and limit the dimension of the diagonalization to 5 · 10 8 . We have also preformed calculations of the stable isotopes 36,38 Ar which are of the same quality as those for 40 Ar. Nonetheless we have not included them in our study, because the natural abundance of these isotopes is very minor with less than 0.3%. Finally, there are five stable stable germanium isotopes 70,72,73,74,76 Ge. Consistently with our study of 73 Ge in Ref. [47] we use the RG effective interaction [88] in a configuration space consisting on the 1p 3/2 , 0f 5/2 , 1p 3/2 , and g 9/2 single-particle orbitals. We have performed calculations with alternative effective interactions in the same configuration space [89][90][91][92], and while the excitation spectra may be somewhat different to the ones predicted by the RG interaction, the impact on the nuclear structure factors is very small at the momentum transfers relevant to direct detection searches. Figures 1-6 compare the low-energy excitation spectra of the stable isotopes of fluorine, silicon, argon, and germanium with our theoretical predictions. In all cases, our calculations are in very good agreement with experiment, especially for nuclei with even number of nucleons. In some cases, especially for the odd-mass nucleus 73 Ge, some experimental states are not well reproduced. This is very likely due to the limitation of the configuration space used in our calculation, because the description of 73 Ge is of similar quality with the other effective interactions we have studied.
To complement the assessment of the quality of the nuclear structure calculations, Tables II and III compare theoretical and experimental electromagnetic observables for all these nuclear targets. The comparison includes nuclear charge radii, electromagnetic moments for ground and lowest-excited states, and nuclear matrix elements for selected electromagnetic transitions between low-lying states. Details on the shell model calculation of  nuclear moments and matrix elements can be found, e.g., in Ref. [75]. Charge radii [97], the properties most relevant for coherent nuclear structure factors, are very well reproduced, better than 3% in all cases, and better than 1% in the heavier argon and germanium. For fluorine and silicon, Table II shows an excellent agreement between theory and experiment, as the majority of the predictions reproduce data within experimental uncertainties. For argon and germanium, Table III also shows a reasonable agreement of the nuclear structure calculations with experiment. Nuclear matrix elements within the same isotope can vary over two orders of magnitude, and the theoretical results reproduce well the corresponding hierarchy. In germanium, some theoretical electric moments and transitions underestimate experiment moderately. This suggests that the configuration space used in the calculation is not sufficient to fully account for the most collective states, consistently with the findings from the comparison of the 73 Ge excitation spectrum. We do not expect, however, any significant effect in the ground states involved in elastic WIMP-nucleus scattering.    [93]. Quadrupole moments Q > 0 (Q < 0) indicate nuclei with prolate (oblate) deformation. Theoretical results are compared to experimental data from Refs. [94][95][96].
In Sec. IV and the Python notebook we also cover nuclear structure factors for xenon, which we studied in the context of coherent SI scattering in Refs. [22,49], and more generally in Refs. [46,47]. The quality of the nuclear structure calculations of stable xenon isotopes is similar to that of argon or germanium, the heaviest nuclear targets considered in this section. Calculations for xenon have been compared to experimental data in Refs. [22,47] for excitation spectra, and in Refs. [46,98] for electromagnetic properties. Given that the charge radii are closely connected to the nuclear structure factors considered here, Table IV compares theoretical results with experiment. In all isotopes the calculations reproduce measured radii to better than 1%.

IV. STRUCTURE FACTORS
Apart from the dependence on q, m N , and m χ that is predicted by ChEFT, the generalized SI WIMP-nucleus cross section in Eq. (1) depends on six independent structure factors F. Four of them correspond to the coupling of the WIMP to one nucleon, which can be the same for protons and neutrons (as in the two isoscalar structure factors) or opposite (as in the two isovector ones). In addition, two independent structure factors characterize the simultaneous coupling of WIMPs to two nucleons. In this section we evaluate these six structure factors for the nuclear targets considered in present and future direct detection experiments. An overview of the various contributions, excluding interference terms, is provided   in Fig. 7. In the following, we present in detail our results for the one-body (1b) and two-body (2b) structure factors.
A. One-body structure factors As discussed in Refs. [34,35] there are two different nuclear responses describing the coupling of a WIMP to a single nucleon, F M and F Φ , which receive coherent contributions from several nucleons in the nucleus. In addition to the standard SI scattering, the nuclear response F M describes a subleading contribution that corresponds to the NREFT operators O 5,8,11 , but by far the most important among them is O 11 . In addition, the so-called radius correction to the standard SI structure factor is also coherent [49]. Dropping the contributions from O 5,8 , the scattering cross section including one-nucleon couplings simplifies to wherec M I can be identified with c M,11 I from Eq. (1). Even though there are only four independent isoscalar contributions (plus four isovector ones), in the most general case where all contributions in Eq. (20) compete, the interference of all of them generates a plethora of individual terms that could be considered. Figures 8, 9, 10, and 11 show the leading contributions to the cross section for 132 Xe, 74 Ge, 40 Ar, and 19 F, respectively. The results for xenon use the results of Ref. [49]. Figures 8, 9, 10, and 11 assume that all couplings c are equal to 1, and for the NREFT O 11 term m χ = 2 GeV, which implies that for heavier WIMPs the importance of this term will always be smaller than in the figures. Figures 8, 9, 10, and 11 highlight that the standard SI contribution (solid black lines) is indeed expected to be dominant. Moreover, leading corrections to generalized SI scattering come from the interference of the standard SI and other terms, such as its isovector counterpart (dotted-dashed black), the isoscalar radius corrections (solid blue), or the isoscalar F Φ term (solid red). The only exceptions are, first, the purely isovector SI structure factor (dashed black), suppressed by one or two orders of magnitude. Second, the O 11 contributions (solid green), suppressed by around four orders of magnitude in all cases (the suppression will be larger for heavy WIMPs m χ > 2 GeV).
The variable importance of these contributions is set by the nuclear structure of the corresponding nuclear targets. Isovector contributions are relatively more important in neutron-rich xenon than in the N ≈ Z fluorine. On the other hand, the F Φ contributions are relatively more important in the heavier xenon and germanium, because these targets have more nucleons in single-particle orbitals with aligned spin and orbital angular momentum. In contrast, the F Φ contributions are more suppressed in lighter targets such as argon and especially fluorine, which tend to have nucleons more equally distributed in orbitals with spin parallel and antiparallel to the orbital angular momentum.

B. Two-body structure factors
The two-body amplitudes of the three scalar channels, scalar-scalar (SS), trace anomaly θ µ µ ("θ" in short), and spin-2 (referred to by "(2)"), read where the couplings f π , f θ π , f (2) π are defined in Eq. (8). σ i and τ i refer to spin and isospin operators for nucleon i, F π is the pion decay constant, and g A the axial charge of the nucleon. Throughout, we use PDG values [99], except for the particle masses, for which we use isospin averages of m N = 938.92 MeV and M π = 138 MeV.
While Refs. [43,45,49] introduced coherent two-body currents, the present work includes for the first time the contact-term contributions to the θ µ µ response involving C S and C T and the entire spin-2 contribution. In particular the consistent inclusion of the contact operators is a crucial improvement in this work (the result for θ µ   was already used in Ref. [50]), see Sec. V for an extended discussion. By including the relativistic corrections of subleading one-body terms, both in the θ µ µ and spin-2 channels, it is possible to write the three physical responses in terms of just two new structure factors: see Sec. V and Ap-  pendix C for more details and the precise definition of contact operators. The two structure factors in the naive (non-interacting) shell model read structure factors 40 Ar structure factors FIG. 11: Structure factors for 19 F, 1b contributions only including interference terms. The description is as in Fig. 8.
where the second structure factor is normalized as F b (0) = −2E b /M π with the binding energy of the nucleus E b < 0. F π (q 2 ) corresponds to the scalar-scalar two-body current, and by defining F b (q 2 ) as above, the three physical channels (SS, θ, (2)) are described in terms of these two structure factors because F θ π (q 2 ) = 2F π (q 2 ) + F b (q 2 ), 13: Structure factors for 74 Ge from contributions including 2b terms only. The description is as in Fig. 12. Figures 12, 13, 14, and 15 show the structure factors for 132 Xe, 74 Ge, 40 Ar, and 19 F that include contributions from the coupling to two nucleons. Scalar couplings are described by the F π contributions (thick solid orange line), which can interfere with the SI contribution structure factors 40 Ar 14: Structure factors for 40 Ar from contributions including 2b terms only. The description is as in Fig. 12. structure factors (thin solid orange) and its isovector counterpart (dotteddashed orange) and with an independent two-nucleon coupling (maroon dotted-double-dashed).
The two-nucleon coupling to the trace anomaly receives two contributions. According to Eq. (25), the first one can be described by F π and the second one by the structure factor F b . Figures 12, 13, 14, and 15 show the F π , F b structure factors, and their interferences with one-nucleon couplings. In particular, besides the terms described above, the figures show the full F b structure factor (thick solid indigo line), its interference with the

SI term (thin solid indigo), and its isovector counterpart (dotted-dashed indigo), and finally the interference with the radius correction term (dotted indigo).
Similarly, spin-2 two-nucleon couplings contribute via F π and F b terms according to Eq. (25). Figures 16 and 17 compare the contribution of the physical combination of the F π and F b structure factors that originate in the scalar, trace anomaly, and spin-2 two-nucleon couplings. Taking all coefficients c to unity, the dominant effect is given by the trace anomaly (dashed line), followed by the spin-2 term (dotted-dashed), and the scalar contribution (solid). However, we stress that this hierarchy only reflects the nuclear structure aspects and does not need to be followed by particular models with definite Wilson coefficients. The nucleon matrix elements and especially the BSM couplings can alter the importance of the different structure factors in Figs. 16 and 17.

V. CONTACT OPERATORS IN TWO-BODY STRUCTURE FACTORS
The results for the two-body currents presented in Sec. IV B are based on the ChEFT formalism developed in Ref. [49], see Appendix C for more details. In particular, the chiral power counting follows the proposal from Ref. [100,101], in which the scaling of operators is estimated by dimensional analysis. In this way, the leading contribution for a scalar current stems from pion-exchange diagrams, since (N † N ) 2 contact operators require an additional insertion of a scalar source that is counted in the same way as the quark mass matrix and therefore only appears at subleading order, see, e.g., Ref. [102]. In alternative formulations [103,104] where contact operators are promoted to lower orders, the pion-exchange diagrams would be accompanied by additional contact-term contributions at the same order, as also suggested by renormalization group arguments for external currents [105]. A conclusive test of the importance of these contact operators would require detailed studies of the scalar current in light nuclei along the lines of Ref. [59], using recent precision chiral potentials [106][107][108], contrasted to nuclear σ-terms from lattice QCD [109,110]. Work along these lines is in progress.
In practice, the question arises as to how to deal with such potential contact operators in nuclear many-body calculations. In fact, even in the Weinberg power counting contact operators do occur at leading order in the coupling to the trace anomaly θ µ µ and a spin-2 source. These terms were neglected in Ref. [49], so that the resulting θ µ µ response was incomplete (their contributions were, however, included in Ref. [50]). Here, we show how in the Weinberg power counting the contact operators in these channels are canonically renormalized in terms of nuclear binding energies, a mechanism that no longer applies once additional contact terms are introduced.
One-loop corrections have been worked out in Ref. [111]. Likewise the diagrams in Fig. 18 for the nucleon give with loop integral In the two-nucleon sector, see Fig. 19, the first term in Eq. (22) follows from the pion-exchange diagram (a) by means of Eq. (26) (diagram (b) only enters at higher orders). An additional contribution arises from the contact-term Lagrangian: where the second term derives from the NR expansion of the spin vector S µ = (0, σ/2). The corresponding term in the trace gives the second term in Eq. (22), represented by diagram (c) in Fig. 19.
In contrast to the scalar current, the pion-exchange piece in Eq. (22) behaves as a contact term for q i → ∞, due to the momentum dependence of the pion coupling, see Eq. (26). As first noted in Ref. [50], for vanishing momentum transfer these two pieces combine to the N N  potential V N N . Together with the kinetic-energy operator T , this suggests the renormalization prescription where E b < 0 is the binding energy of the nucleus, whose wave function |Ψ should be obtained from N N interactions only. In practice, the kinetic-energy operator formally enters at higher orders but arises naturally from relativistic corrections, while 3N forces correspond to higher-order corrections. In order to determine C S and C T we therefore use experimental binding energies, corrected for Coulomb interactions according to Ref. [112]. Details of the renormalization both for θ µ µ and spin-2 terms are given in Appendix C.
For the numerical analysis we ignore C T contributions, because the magnitude of C T is expected to be small due to approximate SU (4) symmetry of N N interactions at low energies, with corrections that can be shown to be suppressed by 1/N 2 c [113][114][115]. The resulting C S values in Table V are consistent with the expectation from dimensional analysis [116]: where we have used that C T = 0 impliesC1 S0 =C3 S1 . Moreover, the values in Table V also agree with typical fits to the N N system, e.g., at LO Ref. [83] finds C S = (−56.5 . . . − 118.3) GeV −2 for cutoffs in the range R = (0.8 . . . 1.2) fm.

B. Spin 2
Despite only entering at dimension-8 in Eq. (3), spin-2 contributions become relevant, for instance, in the context of heavy WIMPs, where significant cancellations with spin-0 terms have been observed [117], enhancing the importance of higher-order corrections. The relevant operators are the traceless parts of the energy-momentum tensor, given in Eq. (6). As can be seen from the WIMP part of L 2,NR in Eq. (23). An extension to subleading components is straightforward, since due to Lorentz invariance the pion-exchange contribution becomes proportional to q µ q ν − gµν 4 q 2 . Then the full expression can be reconstructed from the 00 component, which we have identified with F (2) π (q 2 ). The chiral realizations of the matrix elements ofθ µν q,g have been studied in detail in the literature, both for the pion and the nucleon [118][119][120][121][122][123][124]. Here, we only retain the leading couplings related to moments of pion and nucleon PDFs, resulting in the one-and two-body contributions M 1,NR and M (2) 2,NR , as well as the relativistic corrections in Eq. (C6). Motivated by the EMC effect [125], similar methods have been applied in the context of spin-2 couplings in multi-nucleon systems [126,127]. Therefore, measurements of nuclear PDFs could, in principle, provide independent cross checks on the resulting spin-2 structure factor. However, in practice this is not possible with currently employed parameterizations [128][129][130][131]: in the dark matter context, the main effect of the two-body corrections modifies the normalization at q = 0 away from the fully coherent single-particle expectation F(0) = A, see [45,49] for the scalar channel. Presently, nuclear PDFs q A (x) are studied based on bound-proton PDFs q p/A (x) restricted onto the range x ∈ [0, 1], in such a way that the full PDF is reconstructed by q A (x) = Zq p/A (x) + N q n/A (x). The moments of this q A (x) are therefore, by definition, normalized to the coherent limit A and cannot be used to cross check the twobody effect that we derived from the spin-2 couplings of the pion.

VI. SUMMARY
We have presented a comprehensive analysis of the generalized SI scattering of spin-1/2 and spin-0 WIMPs off atomic nuclei. Our analysis considers all contributions that can receive the coherent enhancement from several nucleons in the nucleus, keeping terms up to third order in ChEFT. This includes both the coupling of WIMPs to one nucleon as well as to two nucleons. For two-body interactions we provide, for the first time, a full and consistent treatment of the contact operators that appear at the same order as pion-exchange diagrams, arguing that these contributions can be renormalized to the nuclear binding energy. As a result, just two nuclear structure factors are enough to characterize two-body interactions via scalar operators, the trace of the energy-momentum tensor, and spin-2 operators.
Taking into account all these contributions, we give all one-body and two-body nuclear structure factors relevant for the coherent WIMP scattering off fluorine, silicon, argon, germanium and xenon, covering the targets of the most advanced direct detection searches. For that purpose, we perform large-scale nuclear shell model calculations with configuration spaces and nuclear interactions that describe very well the structure of these nuclei.
Our analysis identifies the parameters that can, at least in principle, be separately constrained in direct detec-tion experiments. These parameters subsume both the BSM couplings of WIMPs with quarks and gluons and the hadronic matrix elements that embed these quarklevel operators into hadrons. The corresponding matching relations are illustrated in detail for both spin-1/2 and spin-0 WIMPs.
The main results of our work, encoded in the nuclear structure factors and the relation between directdetection experiments and BSM couplings, are available as supplementary material in a Python notebook. These results form the basis for a comprehensive study of WIMP-nucleus interactions based on ChEFT. Future extensions concern non-coherent WIMP-nucleus interactions, for which more parameters and nuclear structure factors need to be considered. Accordingly, if the coherent contributions studied in this paper are strongly suppressed, the identification of the underlying quarklevel interactions becomes even more challenging. On the other hand, progress in ab initio nuclear theory paves the way towards fully consistent structure factors from many-body calculations based on ChEFT [51,59,60,[77][78][79]. Such improved nuclear structure factors, including their momentum-dependence, will further help distinguish among possible BSM scenarios.
For the matching onto NR single-nucleon operators we use the conventions with and as well as the operator basis where q = |q|. The expressions above refer to the WIMP-nucleon system. In the nucleus, the operator v ⊥ generates two kinds of contributions [34,35]. First, there are operators dependent on the WIMP velocity with respect to the nucleus, v ⊥ T [see the terms involving ξ i (q, v ⊥ T ) in Eq. (1)]. These terms are very suppressed in the scattering ampli- On the other hand, v ⊥ generates terms that contain the nucleon velocity operator. In this case, the operators are mildly suppressed by q/m N , and include an additional derivative. The latter terms are fully responsible for the F Φ structure factor.
Back to the WIMP-nucleon level, the coherently enhanced terms listed in Table I, complemented by the leading SD response, are derived from where we have ignored the non-coherent terms in the NREFT expansion. For the remaining couplings the momentum dependence is indicated by the relativistic momentum transfer t, which reduces to t = −q 2 up to relativistic corrections. Some of the amplitudes in Eq. (A5) receive contributions that break Galilean invariance. For such terms, which only appear beyond O(p 3 ) in the chiral expansion, we assume center-of-mass kinematics, for which the velocity in Eq. (A3) simplifies to In principle, corrections to this identification would need to be considered when calculating the nuclear structure factors, similar to the boost correction in Ref. [132], but given that these contributions are already highly suppressed, we only keep the center-of-mass component. The appearance of such Galilean-invariance-breaking terms at subleading orders in the NR expansion has been pointed out in Ref. [58]. However, at variance with Ref. [58], we already find such contributions in the context of the Pauli form factor F 1 in the V V channel. This is reflected by the corresponding coefficients of O 3 and O 5 . We find similar discrepancies to Ref. [58] in the NR expansion of the tensor current. Besides the O 3 and O 5 coefficients, we also disagree in that in our expressions the induced form factors F 2,T and F 3,T [see Eq. (A8) below] combine to the tensor magnetic moments, so that the less well determined individual form factors are not required at this order in the expansion. Expressed in terms of nucleon matrix elements we have q,N + C (2) g f (2) g,N , For the axial-vector and pseudoscalar matrix elements g N A , g N P , and h N 5 we refer to Ref. [43], since in the present paper the numerical analysis is restricted to the coherently enhanced contributions.
The scalar couplings f N q , scalar radiiσ andσ s , as well as charge radii r 2 E N , r 2 E,s N and magnetic moments κ N , κ s N are discussed in detail in Ref. [49] (see also Ref. [133] for the heavy-quark couplings). For strangeness in the vector channel there is increasing evidence from lattice QCD that these couplings are extremely small, we take κ s N = 0.006(4) and r 2 E,s N = 0.0012(9) fm 2 from Ref. [134]. The scalar couplings to u-and d-quarks can be reconstructed [135] from the pion-nucleon σ-term and input for the proton-neutron mass difference [136][137][138][139]. Here, the tension between phenomenology [140][141][142] and lattice QCD [143][144][145][146] already discussed in Ref. [49] persists. More recently, the phenomenological determination from data on pionic atoms [147][148][149][150][151] has been confirmed by an independent extraction from low-energy pion-nucleon cross sections [152], making the resolution of the tension all the more pressing.
The tensor form factors of the nucleon are defined according to [153,154] N (p )|qσ µν q|N (p) where the tensor charges of the proton f T,p q = F q,p 1,T (0) have been recently calculated to high precision in lattice QCD (evaluated at scale µ = 2 GeV) [155]: and the neutron ones follow from isospin symmetry according to The other couplingsκ T,p q = F q,p 2,T (0)+2F q,p 3,T (0), related to the tensor magnetic moments κ T,p q = −2κ T,p q , are less well determined. Lattice QCD information [156] is consistent with an estimate based on analyticity and unitarity of the form factors in analogy to Ref. [157], combining the pion tensor charge [158] with the electromagnetic form factors of the nucleon [159]. The resulting valuesκ T,p u = −1.3(5),κ T,p d = −0.7(3),κ T,p s = 0.00(1) [160] indicate that for u-and d-quarks the induced terms in the tensor decomposition in Eq. (A8) are actually dominant. As in the vector case, the strangeness content is very small. In addition, the known pion tensor charge [158] allows us to calculate the pion-exchange diagram also in this channel, but, similarly to the vector current, this contribution is suppressed by (N − Z)/A due to its isospin structure.
Finally, the spin-2 couplings of the nucleon are given as moments of nucleon PDFs q(x) in complete analogy to Eq. (10): subject to the same sum rule q f (2) q,N + f (2) g,N = 1.

Appendix B: Matching to NREFT for scalar dark matter
For a spin-0, Standard-Model singlet χ the analysis is based on the effective Lagrangian where the notation follows closely the spin-1/2 case in Eq. (3), and for a real scalar we have C V V q = C T T q = 0. We have not included an axial term of the form χ † i∂ µ χqγ µ γ 5 q, because such a contribution reduces to a combination of O 7 and O 10 NREFT operators, and is therefore even further suppressed than the standard SD interaction in the spin-1/2 case. Similarly, we have neglected a tensor operator withqσ µν iγ 5 q.
Most single-nucleon amplitudes come out as in the spin-1/2 case, up to an additional factor of m χ whenever a derivative is required in the effective operator, and the corresponding factors of Λ. To have the reduced amplitudes in the same conventions as for spin-1/2, a factor m χ needs to be removed which leads to where we have kept the same notation for the couplings as in Eq. (A7), with the understanding that the Wilson coefficients therein now refer to the operators defined in Eq. (B1).
In total, the analog of Eq. (7) for spin-0 becomes where now ζ = 1(2) for a complex (real) scalar. At this order in the ChEFT expansion we do not find contributions from O 5,8,11 .
The two-body θ µ µ and spin-2 amplitudes given in Eqs. (22) and (23) can be rewritten according to where we have used that q = −q 1 − q 2 and thus In both cases most of the first term excluding the q 2dependent piece can be related to M SS 2,NR , leading to the F π contribution in Eq. (25). The q 2 -dependent part can be absorbed into a redefinition of F b in Eq. (24) to avoid the introduction of another structure factor. This redefinition does not change the normalization of F b (0) given by the nuclear binding energy. The remaining term involving pion propagators can be expressed in terms of On the other hand, the relativistic corrections to the one-body θ term and the spin-2 channel read ∆M 1,NR = f (2) N 16m 2 N (p 2 2 + p 2 2 − 3p 2 1 − 3p 2 1 )δ(p 1 − p 1 ) In the limit q → 0 these amplitudes become proportional to the kinetic-energy operator T : Summarizing all 1b and 2b contributions we obtain 1+2,NR =f (2) 1,NR .
In the limit q → 0 we have q 1 = −q 2 , so that the momentum transfers in M (C10) Together with the contact terms in Eqs. (C8) and (C9) we recover twice the complete leading-order chiral N N potential V LO N N . In the limit q → 0 both ∆M θ 1,NR and ∆M 1,NR become proportional to T , leading in both cases to the linear combination T + V LO N N , that can be renormalized to the nuclear binding energy at LO order in the chiral expansion. For the θ µ µ current in Eqs. (C5) and (22) this conclusion follows directly from Eq. (C8). On the other hand, the spin-2 one-body contribution in Eq. (C6) carries a coupling f (2) N that, in general, may differ from f (2) π , the coupling of the two-body term in Eq. (23), see Eqs. (A7) and (8). However, due to the sum rules in Eqs. (11) and (A12), when all Wilson coefficients are set equal we have f (2) N /m N = f (2) π /M π , so that these two couplings cancel in the last term of Eq. (C9) to ensure renormalizability at this order in ChEFT. Since the comparison of Eqs. (15) and (A13) shows that the individual spin-2 couplings do not differ much (especially when considering an isospin average for the nucleon), we assume that for spin-2 the kinetic-energy operator aligns as required. Similarly, in Eq. (23) we also assume that the coefficients of the spin-2 (N † N ) 2 contact operators can be determined in the same way, at least up to higherorder effects.
Structure factor F Si Ar Ge Xe  In practice, once the contribution at q = 0 is renormalized to the nuclear binding energy at LO by adjusting the contact terms C S and C T , the contributions from ∆M θ 1,NR and ∆M (2) 1,NR to the full structure factors are small. In addition, while in principle the structure factors could differ for finite q, due to the different functional form of Eqs. (C5) and (C6) for q = 0, we find that the results are practically the same using either expression. In view of the uncertainties from higher orders in the chiral expansion, this shows that the definition of a single new structure factor is sufficient. We choose N N + 2(C S + C T σ 1 · σ 2 ) which in the limit q → 0 reduces to The corresponding response functions in the naive shell model are F π (q 2 ) = 1 2 occ N 1 N 2 |(1 − P 12 )| 1 f π M SS 2,NR |N 1 N 2 , where due to Eq. (C12) we haveF b (0) = −2E b /M π with the binding energy of the nucleus E b < 0. In Eq. (24) we redefine F b (q 2 ) =F b (q 2 ) − q 2 M 2 π F π (q 2 ). We evaluate the two-body structure factors in Eq. (C13) in a non-interacting shell model, as described in detail in Ref. [49], with the multidimensional integrations performed using the CUBA library [163]. For orbitals in the configuration space of the nuclear shellmodel calculations, the occupation numbers are taken from the full diagonalization discussed in Sec. III. The that depend on the mass number A of the nucleus. The functional form in Eq. (C14) proves an efficient representation of the two-body structure factors as well, but the optimal maximal power needs to be determined empirically. Table VI provides the complete list. second part of the notebook, users can set specific values for the Wilson coefficients that describe the WIMPquark/gluon couplings. The notebook generates the corresponding nucleon and pion matrix elements. Finally, the package yields the response including all channels that contribute to the choice of Wilson coefficients. For completeness, we also implemented a routine that calculates the rate corresponding to the standard halo model (SHM) [21]. In this way we clarify conventions and facilitate the use of improved astrophysical input in the future, as is increasingly becoming available with the Gaia mission [164], see, e.g., Refs. [165][166][167]. The SHM is defined by ρ = 0.3 GeV/cm 3 and where z = v esc /v 0 , v esc = 544 km/s, v 0 = 220 km/s, and the Earth's velocity v E = 232 km/s drops out in the normalization. For the operators considered in the notebook one needs where x min = v min /v 0 , η = v E /v 0 , and z = min(z, x min + η), z = z + η − x min . (D5)