From quarks to nucleons in dark matter direct detection

We provide expressions for the nonperturbative matching of the effective field theory describing dark matter interactions with quarks and gluons to the effective theory of nonrelativistic dark matter interacting with nonrelativistic nucleons. We give the leading and subleading order expressions in chiral counting. In general, a single partonic operator already matches onto several nonrelativistic operators at leading order in chiral counting. Thus, keeping only one operator at the time in the nonrelativistic effective theory does not properly describe the scattering in direct detection. Moreover, the matching of the axial--axial partonic level operator, as well as the matching of the operators coupling DM to the QCD anomaly term, naively include momentum suppressed terms. However, these are still of leading chiral order due to pion poles and can be numerically important. We illustrate the impact of these effects with several examples.

different direct detection experiments [19]. The maximal momentum exchange between DM and the nucleus is q max 200 MeV, see Fig. 1. This means that one is able to use chiral counting, with an expansion parameter q/Λ ChEFT 0.3 to organize different contributions in the nucleon EFT for each of the operators coupling DM to quarks and gluons [1,16,[20][21][22][23][24].
In this paper we rewrite the leading-order (LO) results in the chiral expansion of Ref. [1] in terms of single-nucleon form factors. We also extend these results to higher orders in the (q/Λ ChEFT ) 2 expansion up to the order where two-nucleon currents are expected to become important (for numerical estimates of two-nucleon currents see [16-18, 25, 26]).
The starting point is the interaction Lagrangian between DM and the SM quarks, gluons, and photon, given by a sum of higher dimension operators, Here the C (d) a are dimensionless Wilson coefficients, while Λ can be identified with the mass of the mediators between DM and the SM (for couplings of order unity). The sums run over the dimensions of the operators, d = 5, 6, 7, and the operator labels, a. We keep all the operators of dimensions five and six, and all the operators of dimension seven that couple DM to gluons. Among the dimension-seven operators that couple DM to quarks we exclude from the discussion the operators that are additionally suppressed by derivatives but have otherwise the same chiral structure as the dimension-six operators (for the treatment of these operators see [27]).
Here G a µν is the QCD field strength tensor, while G µν = 1 2 ε µνρσ G ρσ is its dual, and a = 1, . . . , 8 are the adjoint color indices. Moreover, q = u, d, s denote the light quarks (we limit ourselves to flavor conserving operators). Note that we include two more dimension-seven operators than in [1], so that we have all the operators included in [30]. The remaining dimension-seven operators coupling DM to quarks are listed in [27], while the effect of dimension-seven operators coupling DM to photons is discussed in [31]. There are also the leptonic equivalents of the operators Q (6) 1,q , . . . , Q (6) 4,q , and Q (7) 5,q , . . . , Q 10,q , with q → . The aim of this paper is to provide compact expressions for the non-perturbative matching at µ 2 GeV between the EFT with three quark flavors, given by Eq. (1), and the theory of DM interacting with nonrelativistic nucleons, given by For each operator the matching is done using the heavy baryon chiral perturbation theory expansion [32] up to the order for which the scattering amplitudes are still parametrically dominated by single-nucleon currents. The relevant Galilean-invariant operators with at most two derivatives are and in addition where N = p, n. At next-to-leading order (NLO) we also need one operator with three derivatives, Our definition of momentum exchange differs from [6] by a minus sign, cf. Fig. 2, while the operators coincide with those defined in [6]. Each insertion of q is accompanied by a factor of 1/m N , so that all of the above operators have the same dimensionality.
This paper is organized as follows: in Section II we give the matching conditions for

II. FERMIONIC DARK MATTER
The hadronization of operators Q 1,q , . . . , Q 10,q , in Eqs. (3)- (9), leads at LO in the chiral expansion only to single-nucleon currents [1]. The resulting matrix elements scale as q ν LO , where [1,27] [Q counting m q ∼ m 2 π ∼ q 2 , and not displaying a common scaling factor. The LO contributions are either due to scattering of DM on a single nucleon (the first diagram in Fig. 3), or on a pion that attaches to the nucleon (the second diagram), or both. The contributions from DM scattering on two-nucleon currents arise at O(q ν LO +1 ) for O (6) 2,q , O 5,q , and O 6,q , at O(q ν LO +2 ) for O (6) 1,q , and at O(q ν LO +3 ) for all the other operators. Up to these orders, the hadronization of the operators Q 1,q , . . . , Q 10,q can thus be described by using form factors for single-nucleon currents.
The form factors are given by . . .
where we have suppressed the dependence of nucleon states on their momenta, i.e. N | ≡ N (k 2 )|, |N ≡ |N (k 1 ) , and similarly,ū N ≡ u N (k 2 ), u N ≡ u N (k 1 ). The momentum ex- The form factors F i are functions of q 2 only.
The axial current, the pseudoscalar current, and the CP-odd gluonic current receive contributions from light pseudoscalar meson exchanges corresponding to the second diagram in Fig. 3. For small momenta exchanges, q ∼ m π , one can expand the form factors in q 2 , as where we kept both the pion and eta poles and denoted the order of the various terms in chiral counting. The coefficients a i , b i , c i are momentum-independent constants. Note that the pion and eta poles for the GG operator are suppressed by q 2 and are thus of the same chiral order as the constant term, b Ñ G . All the other form factors do not have a light pseudoscalar pole and can be Taylor expanded 1 around q 2 = 0, where the prime on F denotes a derivative with respect to q 2 . The values of F The size of the form factors that do not have light-meson poles are, at zero recoil, (only here and in the remainder of the subsection we use the abbreviation q = u, d).
A. Leading order expressions We first give the expressions for the nonrelativistic EFT Lagrangian (10) at LO in chiral counting. In this case we only need the values of a π i , a η i , b Ñ G , and F i (0). In addition to taking 1 We assume that the NLO terms involving chiral logarithms of the form (m 2 π − q 2 ) log(m 2 π − q 2 ) were also expanded in q 2 . This may give an effective expansion parameter q 2 /(Λ EFT ) 2 with Λ EFT between m π and 4πf ; however, numerically it is found to be closer to the latter, see Appendix A.
the hadronic matrix elements of the quark and gluon currents we also take the nonrelativistic limit of both the DM currents and the nucleon currents. The expressions for this last step are collected in Appendix B. The chirally leading hadronization of the dimension-five operators is thus given by with F N 1 (0) = δ pN the nucleon charge, and µ N the nucleon magnetic moment (see also Appendix A 1). The dimension-six operators hadronize as while the hadronization of the gluonic dimension-seven operators is given by The hadronization of the dimension-seven operators with quark scalar currents results in Q and for the tensor operators Q The nonrelativistic operators have been defined in (11)- (17). In the above expressions all the form factors F q/N i are evaluated at q 2 = 0, apart from F q/N P,P and F Ñ G , where one needs to keep the two meson-pole terms in (29) and the first three terms in (30). The corresponding values of coefficients c N i in the nonrelativistic Larangian, Eq. (10), are given in Appendix E. Several comments are in order. First of all, in several cases a single operator describing the DM interactions with quarks and gluons matches onto more than one nonrelativistic operator in Eqs. (11)- (19) already at leading chiral order. This occurs for 10,q = m q (χiσ µν γ 5 χ)(qσ µν q) ∼ where we only show the approximate dependence on the nonperturbative coefficients (here Q p = 1 is the proton charge, while the values of the axial charge ∆q N , the form factors 3 or Q 7,q , and Q  While it is true that the spin-dependent operator O N 4 can arise from the tensor-tensor operator Q 9,q , this contribution would be of two-loop order in a perturbative UV theory of DM. The axial-axial operator Q (6) 4,q , on the other hand, also leads to spin-dependent scattering and will arise at tree level. Therefore it will, if generated, typically dominate over Q (7) 9,q . The spin-dependent scattering it induces arises from both the O N 4 = S χ · S N and O N 6 = S χ · q S N · q operators. While the latter is O(q 2 ) suppressed, it is enhanced by 1/(m 2 π + q 2 ) at the same time, so that in general the two contributions are of similar size (for scattering on heavy nuclei). In this case, again, one cannot perform the direct detection analysis with The same is true for the operators Q 2,q , Q 3,q , and Q 10,q that each match, at leading order in chiral counting, to at least two nonrelativistic operators. A correct LO description of the DM scattering rate cannot be achieved by using only one or the others.
We explore this quantitatively in Section IV, distinguishing also the cases of light and heavy nuclei.

B. Subleading corrections
We discuss next the NLO corrections to the nonrelativistic reduction of the operators (3)- 5,q , O 6,q , for which the two-nucleon corrections arise at O(q ν LO +1 ), and the operator O Starting at subleading order there are terms that break Galilean invariance. This a consequence of the fact that the underlying theory is Lorentz and not Galilean invariant [33].
These corrections involve the average velocity of the nucleon before and after the scattering event, v a = ( k 1 + k 2 )/(2m N ), and lead to ten new nonrelativistic operators listed in Eqs.
The operators that appear at subleading order in the nonrelativistic reduction can have a qualitatively different structure from the ones that arise at LO. For instance, the vectorvector current operator Q (6) 1,q = (χγ µ χ)(qγ µ q) reduces at NLO to At LO one thus has the number operator O N 1 = 1 χ 1 N and no spin dependence, while the expansion to the subleading order gives in addition velocity-suppressed couplings to spin There is only one operator, where this occurs, though. The tensor-tensor operator, Q 9,q = m q (χσ µν χ)(qσ µν q), leads at LO in the chiral expansion to the spin-spin interaction, O N 4 = S χ · S N . At NLO, on the other hand, one also obtains a contribution of the form ∼ q 2 1 χ 1 N , where we do not display the other q 2 -suppressed terms. For heavy nuclei the coherently enhanced contribution from O N 1 scales as A q 2 /(m N m χ ) ∼ O(1) and thus the formally subleading contribution could, in principle, become important in nuclear scattering. Inspection of this particular case, however, shows that there is a relative numerical factor of 16 enhancing the leading contribution. Furthermore the coherent O(q 2 ) term is suppressed by 1/m N m χ and not simply by 1/m 2 N , further reducing its importance for heavy DM masses. As a result the O(q 2 ) terms are numerically unimportant also for the tensor-tensor operator. In contrast, such coherent scattering is important in µ → e conversion, where the m χ supression gets replaced by m µ [34].
A potential concern is that something similar, but with a less favorable result for the numerical factors, could happen for some other operator due to the uncalculated contributions from the nonrelativistic expansion to even higher orders. However, one can easily convince oneself that this is not the case by using the parity properties of quark and DM bilinears. All the relativistic operators in Eq. (1) that are composed from parity-odd bilinears necessarily involve the parity-odd spin operators for single-nucleon currents at each order in the chiral expansion, because one cannot form a parity-odd quantity from just two momenta -the incoming and the outgoing momentum (cf. (B12)-(B17)). Such operators thus never lead to coherent scattering (the argument above may need to be revisited for two-nucleon currents). This leaves us with the operators composed from parity-even bilinears only. Scalar-scalar operators and vector-vector operators lead to coherent scattering already at LO, giving tensor-tensor operator as the only left over possibility. The reduction of the tensor bilinear, Eq. (B16), gives at LO ∼ µναβ v α S β , while at NLO one also gets, among others, the combination v [µ q ν] . The latter does not involve spin and leads to coherent scattering. However, due to numerical prefactors, the latter contribution is still subleading, as was shown above.

III. SCALAR DARK MATTER
The above results are easily extended to the case of scalar DM. 2 For relativistic scalar DM, denoted by ϕ, the effective interactions with the SM start at dimension six, where ellipses denote higher dimension operators. The dimension-six operators that couple DM to quarks and gluons are while the coupling to photons are Here 4,q , with q → . At LO in chiral counting the operators coupling DM to quark and gluon currents hadronize as The expressions valid to NLO in chiral counting are given in Appendix D.
There are a number of qualitative differences between the cases of fermionic and scalar DM. For instance, since scalar DM does not carry a spin there is a much smaller set of operators that are generated in the nonrelativistic limit. This greatly simplifies the analysis. in. At the end of this section, we will also comment on the NLO corrections. The rate R, i.e., the expected number of events per detector mass per unit of time, is given by where E R is the recoil energy of the nucleus, m A is the mass of the nucleus, and ρ χ is the local DM density. The integral is over the DM velocity v in the Earth's frame with a lower bound given by The DM-nucleus scattering cross section dσ/dE R in Eq. (74) is given by The nonrelativistic matrix element squared is [6] 1 2J χ + 1 where J χ = 1/2 is the spin of DM in our examples and J A is the spin of the target nucleus.
The nuclear response function W i depend on momentum exchange, q ≡ | q |. The spinindependent scattering is encoded in the response function W M which, for instance, arises from the matrix element squared of the nuclear vector current. In the long-wavelength limit, q → 0, W M (0) simply counts the number of nucleons in the nucleus giving coherently The response functions W Σ and W Σ have the same long-wavelength limit and measure the nucleon spin content of the nucleus. W ∆ measures the nucleon angular momentum content of the nucleus, while W ∆Σ is the interference term.
These functions roughly scale as where the actual size depends on the particular nucleus and can differ significantly from one nucleus to another. The prefactors R i encode the dependence on the c N i (q 2 ) coefficients, Eq. (10), and on kinematical factors. For instance, the coefficient of the coherently enhanced term is where The sum in Eq. (76) is over isospin values τ = 0, 1 which are related to the proton and neutron coefficients by c 0 can be found in [6]. Using these expressions for R τ τ M together with our expressions for the hadronization of the EFT operators, Eqs. (39)- (54), which give the coefficients c τ i (see Appendix E), we are now in a position to obtain the rates in a DM direct detection experiment assuming a particular interaction of DM with the visible sector.
In the following, when we calculate the scattering rate and plot the bound on the squared UV Wilson coefficients, we restrict the integral over the recoil energy. To approximate the LUX sensitivity region we integrate over E R ∈ [3, 50] keV for Xenon [29]. To approximate PICO's [28] sensitivity we integrated over E R > 3.3 keV for Fluorine -see Figs. 1 and 5.
To obtain total rates for scattering on Xenon, we assume an exposure of 5000 kg·yr which is representative of the next generation two-phase liquid Xenon detectors. Since Xenon has eight naturally occurring stable isotopes, we sum over them weighted by their natural abundances.
Xenon; C the momentum dependence of c τ τ i are expected to be smaller than those of the momentum dependence of W τ τ i . However, because the pseudoscalar hadronic currents contain pion poles, the corrections due to non-zero momentum in the corresponding c τ τ i are of O( q 2 /m 2 π ) and can be large.
The effect of such contributions for scattering on Xenon is shown in Fig. 4. The chirally leading hadronization of the axial-axial operator (χγ µ γ 5 χ)(qγ µ γ 5 q) contains two nonrela- . The latter is momentum suppressed but comes with a pion-pole enhanced coefficient, see Eq. This is compared with the extraction of the bound on C 4,u /2 and a DM mass m χ = 100 GeV. Independent of the DM mass, however, the pion pole is completely absent for C where ∆q i is the axial charge of quark q i and the δm coefficient is the size of isospin breaking for pion exchange and the SU (3)-flavor breaking for eta meson exchange, see Eq. (A42).
Note that isospin breaking is O(1) for the matrix element of the QCD theta term, unlike in any other observable. The importance of isospin-breaking, but pion-pole enhanced contributions is reflected in the DM scattering rates. The bounds on the Wilson coefficients C 3,4 in Fig. 6, obtained with the correct full form factor dependence, are depicted with solid black lines. For weak-scale DM masses they can be even up to an order of magnitude stronger than the bounds obtained by only using the zero recoil form factor, FG(0) (dashed blue lines). Ignoring the leading q 2 -dependence in FG also leads to a large distortion of the shape in dR/dE R as shown in Fig. 5 (right) for the Q 4 operator and m χ = 200 GeV. In this case, there is a visible change in the shape of the differential rate even for scattering on Fluorine, despite small momenta exchanges. The effect is striking for the scattering on Xenon where momenta exchanges are typically larger. For the Q 3 operator, the distortion is slightly smaller, but otherwise comparable to the one shown.
For the Q 4,q and Q (7) 4 operators discussed above and shown for scattering on Xenon in Figs. 4 and 6 respectively, the q 2 dependence in the meson poles is negligible for scattering on Fluorine. To understand this it is useful to consider the differential scattering rate as a function of the recoil energy. This is shown in Fig. 5 for a fixed DM mass of 200 GeV.
For both interactions, the E R spectra for Fluorine do not differ significantly when the q 2 dependence in the meson poles is neglected since a given value of E R results in momentum transfer q 2 /m that is smaller by an order of magnitude in Fluorine than in Xenon.
A qualitatively different example is given in Fig. 7 which shows the bounds on the Wilson coefficient C 3,q as function of m χ for scattering on Xenon and Fluorine. The vector-axial operator, Q 3,q = (χγ µ χ)(qγ µ γ 5 q), Eq. (4), matches onto two non-relativistic operators, O N 7 ∝ S N · v ⊥ and O N 9 ∝ S χ · ( q × S N ). At leading order in chiral power counting, the hadronization of the axial quark-current in Q the DM mass (i.e, two powers in the rate) and thus becomes subleading for larger DM masses.
Since the contributions are correlated yet scale differently with m χ , it is crucial to consider both non-relativistic operators when setting bounds from direct detection experiments (see, e.g., [36]).
The non-trivial interplay between different non-relativistic operators can also be seen in the case of dipole interaction, Q 1 , shown in Fig. 8. This operator matches onto four (39). Out of these, two are coherently enhanced, One expects these two to dominate for heavier nuclei, as shown explicitly for Xenon in Fig. 8 (left). The O N 5 operator is enhanced by an explicit photon pole prefactor, 1/ q 2 , which overcomes the velocity suppression and leads to its dominance over all other contributions. The contribution from the O N 1 operator, on the other hand, is local and is suppressed for heavy DM by a 1/m χ factor. Its contribution is, therefore, relevant only for light DM.
For DM scattering on lighter nuclei, the situation is more involved. The coherent enhance-  9,q where two coherently enhanced operators, O N 1 and O N 5 , appear at NLO in the expansion, while at LO no coherently enhanced operators are present. However, even for Xenon, the coherent enhancement is not enough to compensate for the q 2 /m N m χ suppression accompanied by a relative factor of 1/16, and thus the resulting correction is of O(5%). A similar coherently enhanced contribution appears for Q (7) 10,q operator at O(q 4 ) and is thus completely negligible.

V. CONCLUSIONS
In this article we derived the expressions for the matching of an EFT for DM interacting with quarks and gluons, described by the effective Lagrangian L χ in Eq. (1), to an EFT We worked to next-to-leading order in the chiral expansion, but also discussed separately the expressions for the leading-order matching. At LO the scattering of DM on nucleons only depends on the DM spin S χ , the nucleon spin S N , the momentum exchange q, and the averaged relative velocity between DM and nucleon before and after scattering, v ⊥ . Numerically the NLO corrections are always small, at the level of O( q 2 /m 2 N ) or a few percent, unless one fine tunes the cancellation of LO expressions. This result is nontrivial for the partonic tensor-tensor operator Q 9,q = m q (χσ µν χ)(qσ µν q), since in that case the LO term is spin-dependent, while the NLO corrections contain a spin-independent contribution that is coherently enhanced. In principle this could compete with the LO term. However, due to fortuitous numerical factors, it remains subleading.
While our results were obtained by assuming that the mediators between the DM and the visible sector are heavy, with masses above several hundred MeV, the formalism can be easily changed to accommodate lighter mediators. In this case the mediators cannot be integrated out, but lead to an additional momentum dependence of the coefficients in the nonrelativistic Lagrangian L NR , Eq. (10), and potentially to a modified counting of chirally leading and subleading terms. The details of the latter would depend on the specifics of the underlying DM theory.
As a side-result, our expressions show that from the particle physics point of view it is more natural to interpret the results of direct detection experiments in terms of an EFT where DM interacts with quarks and gluons, Eq. (1). The reason is that several of the partonic operators in L χ match to more than one nonrelativistic operator already at leading order in chiral counting. In such cases it is then hard to justify singling out just one nonrelativistic operator in the analysis of direct detection experimental results.
The situation becomes even more complicated if the partonic operator matches onto several nuclear operators with different momentum dependence, since in the experiments one integrates over a range of momenta. A cautionary example of wider phenomenological interest is the case of the axial-axial partonic operator, Q 4,q = (χγ µ γ 5 χ)(qγ µ γ 5 q), which induces spin-dependent scattering. At leading chiral order this is described by a combination of the O N 4 = S χ · S N and O N 6 ∼ S χ · q S N · q nonrelativistic operators. Naively the latter is momentum suppressed. We find that this is true for DM scattering on light nuclei, such are never generated as leading operators when starting from a UV theory of DM. They enter only as subleading corrections in the scattering rates, and can always be neglected (as can the other nine nonrelativistic operators listed in Appendix C that have already never been considered).
In conclusion, we advocate the use of partonic level EFT basis Eqs. (2)-(9) as a phenomenologically consistent way of interpreting direct detection data. Including all the variations due to quark flavor assignments there are 34 operators in total, which is not much more than the 28 nonrelativistic operators used at present. Moreover, using the partonic level EFT also has the added benefit of providing a simple connection with the use of EFT in collider searches for dark matter, via straight-forward renormalization-group evolution.
Using the strange magnetic moment [37] (see also [38]) one gets, using isospin symmetry, For the slope of F q/N 1 (q 2 ) at q 2 = 0 one obtains [8] F u/p 1 Above we used the definitions for the proton and neutron matrix elements of the electromagnetic current, where J µ em = 2ūγ µ u −dγ µ d −sγ µ s)/3. The Sachs electric and magnetic form factors are related to the Dirac and Pauli form factors, F N 1 and F N 2 , through [40] (see also, e.g., [41]) At zero recoil one has for the electric form factor, G p E (0) = 1, G n E (0) = 0, while the magnetic form factor at zero recoil gives [35], i.e., the proton and neutron magnetic moments in units of nuclear magnetonsμ N = e/(2m N ).
The anomalous magnetic moments are F p 2 (0) = a p , F n 2 (0) = a n . The charge radii of the proton and neutron are defined through (A12)

Axial vector current
The matrix element of the axial-vector current (23) As always, the coefficients for the neutrons are obtained through a replacement p → n, u ↔ d (no change is implied for g A ). We work in the isospin limit, so that ∆u ≡ ∆u p = ∆d n , ∆d ≡ ∆d p = ∆u n , ∆s ≡ ∆s p = ∆s n .
The derivative of the axial form factor at zero recoil is well known for the u − d current.
Using the dipole ansatz [49] gives F A (0)/F A (0) = 2/m 2 A , with m A the appropriate dipole mass. A global average over experimental [50,51] and lattice [46,52] gives for the u −  (14)GeV [46,53] and for the strange-quark current, m s A = 0.82 (21) GeV [46]. This gives or in terms of normalized derivatives while for the strange quark At NLO F q/N P (q 2 ) needs to be expanded to At NLO the residua of the poles change by corrections of O m 2 π,η /(4πf 2 π ) 2 ≈ 0.01 − 0.05. For instance, for the u − d current one has at NLO in HBChPT [54], whereB 2 ≈ −1.0 ± 0.5 is the HBChPT low energy constant, while r 2 to the values in (A26). Note that these are a small correction to the LO expression when the pion pole is present, but can be important when this is not the case.

Scalar current
The scalar form factors F q/N S , Eq. (24), evaluated at q 2 = 0 are conventionally referred to as nuclear sigma terms, where σ N qū N u N = N |m qq q|N , |N and N | represent the nucleon states at rest. Another common notation is σ N q = m N f N T q . Taking the naive average of the most recent lattice QCD determinations [55][56][57], we find The matrix elements of the u and d quarks are related to the pion-nucleon sigma term, defined as σ πN = N |m(ūu +dd)|N , wherem = (m u + m d )/2. A Heavy Baryon Chiral Perturbation Theory analysis of the πN scattering data gives σ πN = 59 (7) MeV [58], in agreement with σ πN = 52(3)(8) MeV obtained from a fit to world lattice N f = 2 + 1 QCD data [59]. Including, however, both ∆(1232) and finite spacing in the fit shifts the central value to σ πN = 44 MeV. We thus use a conservative estimate σ πN = (50 ± 15) MeV. Using the expressions in [60] this gives For corrections of higher order in chiral counting one would also need F q/N S (0). However, these are always of higher order than the two-nucleon contributions, not captured in our expressions.

Pseudoscalar current
In the LO expressions we only need the light meson pole parts of the pseudoscalar form factor, Eq. (25), The residua of the poles are given by where the values of the axial-vector elements, ∆q, are given in (A18) and (A19). Moreover, B 0 is a ChPT constant related to the quark condensate given, up to corrections of O(m q ), by qq −f 2 B 0 . Using quark condensate from [61] and the LO relation f = f π , with f π the pion decay constant, one has B 0 = 2.666 (57) GeV, evaluated at the scale µ = 2 GeV.
In practice, B 0 never appears by itself, but rather as the product B 0 m q which can be expressed in terms of the pion mass and quark mass ratios, The numerical values are obtained using the ratios of quark masses, m u /m d = 0.46 ± 0.05, m s /m = 27.5 ± 0.3 (see the quark mass review in [35]), and the charged-pion mass m π .
At NLO in the chiral expansion, the above expressions for a q/p P,π and a q/p P,η get corrections of O(m 2 π,η /(4πf π ) 2 ). In addition one needs to keep the constant term in the q 2 expansion of the form factor In our numerical analysis we estimate the size of these higher-order corrections by using the while keeping a q/p P,π , a q/p P,η at their LO values. This treatment of NLO corrections is only approximate, but suffices for the present precision. Furthermore, it can be improved in the future.

CP-even gluonic current
The matrix element of the CP-even gluonic current (26) is parametrized by a single form factor F N G (q 2 ). The LO expressions in chiral counting require only its value at zero momentum transfer, The nonperturbative coefficient m G is the gluonic contribution to the nucleon mass in the isospin limit, The trace of the stress-energy tensor, θ µ µ = −9α s /(8π)G µν G µν + u,d,s m qq q, yields the relation where in the last equality we used the values for σ q in (A28) and (A29). While the isospin violation in the σ N q values is of O(10%), this translates to a very small isospin violation in m G , of less than 1 MeV. The value of m G in (A38) thus applies to both N = p and N = n.
For the derivative of F G at zero recoil we use the naive dimensional analysis estimate (A39)

CP-odd gluonic current
The matrix element of the CP-odd gluonic current (27) is related to the matrix elements of the axial and pseudoscalar currents through the QCD chiral anomaly. Namely, a chiral rotation of the quark fields, q → exp(iβγ 5 )q, shifts the QCD theta spurion by θ → θ −2 Tr β, along with corresponding changes in the pseudoscalar and axial-vector spurions (see Ref. [1]).
This implies a relation, valid at leading order in the chiral expansion. To shorten the notation we defined 1/m = (1/m u + 1/m d + 1/m s ). In terms of form factors this gives The leading order contributions from F q/N P cancel in the sum, giving The pion pole contribution would vanish in the exact isospin limit. However, the isospin breaking effects in the matrix element ofGG operator are not small [62]. This is unlike The LO expression for F Ñ G , Eq. (A42), contains both the constant term as well as poles of the form ∼ q 2 /(m 2 π − q 2 ). At NLO in chiral counting one also has in addition the O(q 2 ) contribution, At NLO the a Ñ G,π , a Ñ G,η , b Ñ G coefficients differ from their LO values in (A42) by relative correction of the size O m 2 π,η /(4πf π ) 2 , while the NDA estimate for the NLO coefficient is c Ñ G ≈ 1.

Tensor current
The matrix element of the tensor current (28) is described by three form factors, F q/N T,0 (q 2 ), F q/N T,1 (q 2 ), F q/N T,2 (q 2 ). These are related to the generalized tensor form factors through (see, e.g., [63,64]) The tensor charges are related to the transversity structure functions δq N (x, µ) by g q T (µ) = 1 −1 dxδq N (x, µ). These structure functions can, in principle, be measured in deep inelastic scattering, but this determination is not very precise. Recent lattice calculations include both connected and disconnected contributions and give, in the MS scheme at µ = 2 GeV [65,66], This agrees well with previous, less precise, determinations [63,[67][68][69][70][71][72][73]. It is interesting to compare (A48) with the results from the constituent quark model [74], g u T = 0.97, g d T = −0.24, as we will have to use this model below. In the nonrelativistic quark model, on the other hand, using just SU (6) spin-flavor symmetry, one gets g u T = 4/3, g d T = −1/3, see, e.g., [75].
The zero recoil values of the other two form factors, F q/N T,1 (0) and F q/N T,2 (0), are less well determined. The constituent quark model of [74] gives We assign a 50% error to the above estimates, taking as a guide twice the difference between the determination of g q T in this model and in lattice QCD (A48). For the s quark we use the very rough estimates The linear combination is in fact much better known thanÃ T,10 (0) and B T,10 (0) separately. The tensor magnetic moments, κ q T , for the u and d quarks were determined using lattice QCD to be, at µ = 2 GeV [76], (no uncertainty is given in this reference). In the constituent quark model of [74] one gets κ u T ≈ 2.0, κ d T ≈ 1.2, which agrees with (A53) within the assigned 50% uncertainty (larger values κ u T = 3.60, κ d T = 2.36 are obtained with a simple harmonic oscillator wave function [74,77]). For the strange quark one obtains from the SU(3) chiral quark-soliton model [78] motivating the ranges in (A51) (in [79] a much smaller value κ s T ≈ 0.01 was found. In Refs. [63,72,80], lattice QCD results for the q 2 dependence of F q/N T,0 for u and d quarks were presented. Averaging over them gives where the errors reflect the differences between the three determinations. For the s-quark For the other two form factors an estimate of the derivative at zero recoil can be made using the results from the constituent quark model of [74], giving These estimates most probably have large errors, since within this model one gets T,0 (0) ≈ 0.24 GeV −2 , about a factor of three smaller than lattice QCD determination in (A55). For the strange quark form factor we vary the derivative at zero recoil in the range motivated by the slope dκ s T /dq 2 ≈ −2.2 GeV −2 that one can deduce from the results in [79].

Appendix B: Nonrelativistic expansion of currents for fermions
In this appendix we give the nonrelativistic expansion of the DM and nucleon currents.
We first focus on fermionic DM and then translate the results to nonrelativistic nucleons.
In order to get rid of the time derivative, v · ∂, in the higher-order terms in the Heavy Dark Matter Effective Theory (HDMET) Lagrangian, the tree level relation is supplemented with a field redefinition 4 [82] χ In this way one obtains the conventional "NRQED" Lagrangian, 4 In order for the scattering rates to be independent of this arbitrary field redefinition, contributions to the scattering amplitude from the time-ordered product of the Lagrangians (10) and (B3) have to be included [81]. An explicit calculation shows that, with our choice (B2), these additional contributions vanish to O(p 2 ).
also beyond O(p 2 ) order.
Using (B1) together with (B2) and applying the equation of motion derived from Eq. (B3) we obtain for the DM currents χ v , and S µ = γ µ ⊥ γ 5 /2 is the spin operator. The square brackets in the last line denote antisymmetrization in the enclosed indices, while the ellipses denote higher orders in 1/m χ . We also used the relation where µναβ is the totally antisymmetric Levi-Civita tensor, with 0123 = 1, and The same expressions apply also for nucleon currents, with the obvious replacement χ → N .
In terms of the momenta instead of derivatives the expansions arē where we used the shorthand notation p µ 12 = p µ 1 + p µ 2 . The corresponding expansion of the nucleon currents is obtained through the replacements χ → N , These terms signal that the underlying theory is, in fact, Lorentz rather than Galilean invariant.
In addition to the nonrelativistic operators (11)- (19) there are three new operators of four new operators of O(q 2 ), and three of O(q 3 ), Next we give the expressions for the nonrelativistic reduction of the operators (3)-(9) to subleading order in q 2 . For each of the operators we stop at the order at which one expects the contributions from the two-nucleon currents. We explicitly include a factor in order to convert from the usual relativistic normalization of states, χ(p )|χ(p) = 2E p (2π) 3 δ 3 ( p − p), where E p = p 2 + m 2 χ , to the normalization used in [6]. The hadronization of the dimension-six interaction operators, including the subleading orders for singlenucleon currents, are then given by, Q The terms in the curly brackets arise for the first time at subleading order, i.e., at O(q ν LO +2 ).
The form factors in these expressions are evaluated at q 2 = 0, i.e., F i → F i (0). In the LO terms, on the other hand, one should expand the form factors to O(q 2 ), i.e., in the expressions outside curly brackets, Note that the hadronization of Q 1,q is expected to receive contributions from two-nucleon currents at O(q 2 ), i.e., at the same order as the displayed corrections from the single-nucleon current. In the hadronization of Q (6) 2,q we do not show the subleading corrections from expanding the single-nucleon currents. In this case the two-nucleon currents enter at O(q 2 ), while the higher-order corrections from single-nucleon currents start only at O(q 3 ). Note also that, at O(p 4 ), the hadronization of Q (6) 4,q receives a contribution that is coherently enhanced, but suppressed by a numerical factor ∼ 1/(16m N m χ ).
The hadronizations of the dimension-seven operators are given by 5,q →F q/N S The expressions that appear for the first time at O(q ν LO +2 ) are collected inside the curly brackets. In these the form factors are to be expanded to LO in chiral counting, as denoted in Eqs. (29)- (31). In particular, the form factors without light meson poles are evaluated at q 2 = 0, i.e., for these F i → F i (0) inside curly brackets. In the terms outside curly brackets, however, the form factors should be expanded to NLO, cf. Eqs. (29)- (31). The operators Q 5,q and Q 6,q receive contributions at O(q ν LO +1 ) from two-body currents, so we do not display the corrections from expanding the single-nucleon currents which, in this case, start at O(q ν LO +2 ).

Appendix D: Nonrelativistic expansion for scalar DM
To derive the HDMET for scalar DM, we factor out 5 the large momenta, followed by a field redefinition where the non-relativistic operators are defined in Eqs. (11)- (19) and Eqs. (C2)-(C7).
The coefficients for neutrons are obtained by replacing p → n, u ↔ d. Above we kept only the chirally leading contributions and listed the results only for the non-vanishing c N NR,i (i.e., one has c N NR,2 = c N NR,3 = 0). For the coefficient c N NR,1 , we also kept the q 2 -suppressed contribution fromĈ 9,q , (7) that is, however, coherently enhanced. The contributions due to the magnetic and electric dipole operators, Eqs. (2), are given in Appendix A of [1].