Non-relativistic Effective Interactions of Spin 1 Dark Matter

We investigate the non-relativistic reduction of simplified models for spin 1 dark matter (DM) with the aim of identifying features in the phenomenology of DM-quark interactions which are specific to vector DM. In the case of DM-quark interactions mediated by a spin 1 particle, we find two DM-nucleon interaction operators arising from the non-relativistic reduction of simplified models for spin 1 DM that are specific to spin 1 DM, and which were not considered in previous studies. They are quadratic in the momentum transfer, linear in a symmetric combination of polarisation vectors for the DM particle, and arise from simplified models which do not generate momentum transfer independent operators as leading interactions in the non-relativistic expansion of DM-nucleon scattering amplitudes. Within these simplified models, the new operators cannot be neglected when computing DM signals at direct detection experiments. For example, we find that nuclear recoil energy spectra computed by including or neglecting the new operators can differ by up to one order of magnitude for nuclear recoil energies larger than about 20 keV and DM masses below $50$ GeV. Furthermore, the shape of the expected nuclear recoil spectra depends significantly on whether the new operators are taken into account or not. Finally, neglecting the contribution to DM direct detection signals from the new operators leads to inaccurate conclusions when assessing the compatibility of a future direct detection signal with CMB constraints on the DM relic density, especially when the number of signal events is small, e.g. $\mathcal{O}(1)$.


Introduction
Measurements of velocity dispersion in stellar and galactic systems, gravitational lensing events in galaxy clusters, the hierarchical formation of large scale cosmological structures, and the second peak in the angular power spectrum of the cosmic microwave background radiation show that invisible mass, or Dark Matter (DM), must be present in our Universe [1]. While the microscopic constituents of DM have so far escaped detection, complementary experimental methods are currently used to search for the elusive "DM particle" [2]. Direct detection experiments primarily search for nuclear recoils induced by the non-relativistic scattering of Milky Way DM particles in low-background experiments located deep underground [3,4]. For a review of current direct detection methods, see for example [5].
Predictions for the expected rate of nuclear recoil events at direct detection experiments have been refined significantly in recent years, e.g. [6]. Special emphasis has been placed on modelling the "elementary" DM-nucleon interactions and computing the associated nuclear response functions [7][8][9][10]. Three approaches have been pursued to model the interactions of DM with nucleons: 1) non-relativistic effective theories for DM-nucleon interactions [11,12]; 2) relativistic effective field theories [13]; 3) and, finally, simplified models [14][15][16][17]. In non-relativistic effective theories the relevant degrees of freedom are DM and nucleons. Their interactions are constrained by Galilean invariance and built in terms of a set of basic non-relativistic quantum mechanical operators (see Eq. (3.1)). Relativistic effective field theories for DM direct detection are built in terms of higher dimensional DMquark and -gluon interaction operators constrained by Lorentz invariance. Their matching onto nucleon-level non-relativistic operators has been performed by using chiral effective field theory [18,19]. Not all operators appearing in non-relativistic effective theories for DM-nucleon interaction arise as leading operators through this matching procedure. Finally, simplified models for DM extends the Standard Model of particle physics by one DM candidate and one additional particle that is responsible for the interactions of DM with quarks and gluons. For a discussion on the constraints unitarity and anomaly cancellation impose on this latter approach to DM model building, see for example [20][21][22]. The three approaches are related as follows. When the momentum transferred in DM-nucleus interactions is smaller than the mass of the particle mediator, simplified models reduce to relativistic effective field theories. When the momentum transferred is smaller than the nucleon mass and DM moves at non-relativistic speeds, each relativistic DM-quark and -gluon operator reduces to a linear combination of non-relativistic quantum mechanical operators for DM-nucleon interactions. In the case of spin 0 and spin 1/2 DM, the relation between simplified models and non-relativistic quantum mechanical operators for DM-nucleon interactions has been investigated extensively in recent years [23][24][25]. On the other hand, the non-relativistic reduction of simplified models for spin 1 DM has been significantly less explored.
This work focuses on the non-relativistic reduction of simplified models for spin 1 DM, exploring the cases of a spin 0 and spin 1 mediator particle separately. For spin 1 mediators, we find that two DM-nucleon interactions operators that were not considered in previous studies can actually arise from the non-relativistic reduction of simplified models for spin 1 DM. We also find that the new operators, denoted here by O 19 and O 20 , can have an important impact on the calculation of DM-nucleus scattering rates, especially for DM masses below 50 GeV. In this mass range, scattering rates computed by including O 19 and O 20 or neglecting them differ by up to one order of magnitude in the high energy tail of the nuclear recoil energy spectrum. Focusing on a simplified version of DARWIN [26], we also show the importance of taking into account the new operators O 19 and O 20 when assessing the compatibility of a direct detection signal with the CMB constraint on the DM relic density.
The article is organised as follows. In Sec. 2 we introduce the simplified model framework studied in this work. Sec. 3 focuses on the non-relativistic reduction of the simplified models introduced in Sec. 2 2 Simplified models for spin 1 DM This section introduces the theoretical framework that we will use to investigate nonrelativistic DM-nucleus scattering (Sec. 4) and DM chemical decoupling (Sec. 5) in the case of spin 1 DM. We will consider two cases separately. In the first one, a spin 1 DM candidate interacts with quarks through the exchange of a spin 0 mediator particle, as described below The Lagrangian corresponding to this first case, L φ , extends the Standard Model Lagrangian, L SM , by kinetic and interaction terms involving the vector field X µ describing the spin 1 DM candidate, and the scalar (or pseudoscalar) field φ associated with the spin 0 mediator. In Eq. (2.1), X µν is the DM field strength tensor, and a dagger denotes Hermitian conjugation. Without loss of generality, we assume the coupling constants b 1 , b 2 , h 1 , h 2 , µ 1 , µ 2 and λ X to be real, and denote DM and mediator mass by m X and m φ , respectively. In Eq. (2.1), the spinors q describe the Standard Model quarks. A summation over flavour indexes is understood, and quark interactions are assumed to be flavour diagonal and of universal strength. The 4 × 4 matrices γ µ , with µ = 0, . . . , 3, and γ 5 are the familiar γ-matrices.
In addition to Eq. (2.1), we also consider the case of a spin 1 DM candidate interacting with quarks through the exchange of a spin 1 mediator particle, as described by the following Lagrangian where m G and G µν are, respectively, mass and field strength tensor of the spin 1 mediator G µ . Without loss of generality, we assume the coupling constants λ G , b 3 , b 4 , b 5 , h 3 and h 4 in Eq. (2.2) to be real, while b 6 and b 7 are required to be complex, since their phase cannot be reabsorbed by a field redefinition. As for the case of a spin 0 mediator, we assume flavour diagonal and universal interactions, and a summation over quark flavours is understood. While the Lagrangians in Eqs. (2.2) and (2.1) have already been used in the study of DM particle phenomenology at the LHC [16,17,27], and at direct detection experiments employing ton-scale [28,29] or spin-polarised targets [30], in Sec. 3 we will show that Eqs. (2.2) predicts new, as yet unexplored DM-nucleon interactions in the nonrelativistic limit. As we will see, the new interactions can have an impact not only on the expected rate of DM-nucleus scattering events at direct detection experiments (Sec. 4), but also on the compatibility between a DM direct detection signal and the CMB constraint on the DM relic density (Sec. 5).

New non-relativistic interactions for spin 1 DM
The non-relativistic interaction between DM particles and nuclei can be described in terms of nucleon-level quantum mechanical operators [11,12]. In the most general case, these interaction operators can be constructed from the Hermitian and Galilean-invariant objects [16,17,27]: where 1 N and 1 X are identity operators in nucleon and DM spin space, q is the momentum transferred in the DM-nucleon interaction, and v ⊥ ≡ v + q/(2µ) is the transverse DMnucleon relative velocity obeying v ⊥ · q = 0.
Here v and µ are the DM-nucleon relative velocity and reduced mass, respectively. We denote by m N the nucleon mass. The vectors S N and S X are the nucleon and DM spin operators, while the Galilean-invariant S is a symmetric combination of DM polarisation vectors. An explicit expression for S is given below in Eq. (3.4). Tab. 1 shows the operators defining the effective theory of DM-nucleon interactions [11,12]. We obtain the nucleon-level quantum mechanical operators associated with Eqs. (2.1) and (2.2) as follows. We first integrate out the particle mediators (φ in the case of Eq. (2.1) and G ν in the case of Eq. (2.2)), assuming that the momentum transferred in the scattering is much smaller than the mediator mass. Then we match the resulting quark-level Lagrangians on to relativistic nucleon-level Lagrangians for DM-nucleon interactions by using nucleon form factors from [28]. Finally, we use the nucleon-level Lagrangians obtained in this manner to compute the corresponding amplitudes for DM-nucleon scattering, which we expand at second order in |q|/m N and first order in v ⊥ . In the case of a spin 0 mediator, we find while in the case of a spin 1 mediator described by Eq. (2.2), we obtain the non-relativistic amplitude where 2S r r N = ξ †r σ N ξ r , and σ N is a vector whose components are the three Pauli matrixes acting on the nucleon spinors ξ r , r = 1, 2. Here, we omitted the isospin indexes carried by spinors ξ r and used the following compact notation to simplify the above equations, e.g., h N 3 δ r r = ξ †r h N 3 ξ r ≡ ξ †r [h p 3 (1 + τ 3 )/2 + h n 3 (1 − τ 3 )/2] ξ r , where 1 and τ 3 are the identity and third Pauli matrix acting on the nucleon isospin space. The couplings to protons, h p k , and neutrons, h n k , equal the constants h k , k = 1, . . . , 4, in Eqs. (2.1) and (2.2) times a nucleon form factor that we take from [28] (and references therein). In Eq. (3.3), we only considered terms at most quadratic in q and linear in v ⊥ , and expressed the DM particle spin S s s X and the S s s operator in terms of three-dimensional DM polarisation vectors e si and e s j , In a reference frame where the z-axis is in the direction of the three-dimensional DM particle momentum, the four-dimensional DM polarisation vectors, ε sµ (p) and ε s µ * (p ), Table 1. List of quantum mechanical operators defining the non-relativistic effective theory of DM-nucleon interactions [11,12,28], augmented by including the operators O 19 and O 20 that we found to arise from the non relativistic reduction of simplified models for spin 1 DM with a vector mediator (see Eq. (2.2)). Table 2. Non zero coupling constants appearing in Eq. (3.6). They arise from the non-relativistic reduction of Eq. (2.2). The coupling constants h N k , k = 1, 2, 3, 4, with N = n for neutrons and N = p for protons, equal the coupling constants h k times a nucleon form factor that we take from [28] (and references therein). Notice that the coupling constants in this table are a factor of 2 smaller than the ones in [27,29]. There, we build M NR (defined in Appendix A) following the prescription outlined in [28].
can be written as where the sum runs over the non zero coefficients in Tab. 2, and angle brackets denote expectation values between initial and final DM and nucleon spin and isospin states. The amplitude in Eq. (3.2) can also be written as in Eq. (3.6), but in this case only two operators contribute to the sum and the corresponding coupling constants are given by: 19 and O 20 in Eq. (3.6) were not included in previous classifications, and their phenomenology will be investigated here for the first time. Explicitly, they read as follows They are quadratic in the momentum transfer, linear in the symmetric combination of polarisation vectors for the DM particle, S, and specific to spin 1 DM. As we will see in the next sections, they are not negligible in the analysis of DM direct detection experiments, since they arise from simplified models which do not generate momentum transfer independent operators as leading interactions in the non-relativistic expansion of DM-nucleon scattering amplitudes.

Expected rate of DM-nucleus scattering events
The rate per unit detector mass of DM-nucleus scattering events in a direct detection experiment is is the DM velocity distribution in the detector rest frame, ρ χ is the local DM density, m A is the mass of the target nuclei, E R is the nuclear recoil energy, and dσ/dE R is the differential cross-section for DM-nucleus scattering. For the local DM density and velocity distribution in the detector rest frame we assume, respectively, ρ χ = 0.4 GeV cm −3 and where v e (v esc = 544 km s −1 ) is the Earth (escape) velocity in the galactic rest frame, and the most probable speed, v 0 = 220 km s −1 , equals the circular speed of the local standard of rest. The differential cross-section for DM-nucleus scattering in Eq. (4.1) can be expressed as follows where T f i is the matrix element between initial and final DM-nucleus states of the underlying nuclear potential, and a bar denotes average (sum) over initial (final) DM and nuclear spin configurations. We calculate the matrix element T f i for the cases of a spin 0 and a spin 1 mediator separately. To this end, we first calculate the inverse Fourier transforms of and sum the result over all constituent nucleons to obtain the position-space nuclear potential associated with the DM and mediator model under investigation. Then we evaluate the matrix element between initial and final DMnucleus states of the obtained nuclear potential to obtain T f i . Following [7], for T f i we find where initial and final nuclear states are denoted by |i and |f , respectively. Nuclear charge operators (S τ and T τ ) and nuclear vector operators (P τ , Q τ , and R τ ) are defined in Appendix B. They arise from the expansion in spherical harmonics of the matrix element between initial and final DM states of the position-space nuclear potential obtained from M (0) and M (1) . Summation over the constituent nucleons is encoded in the definition of the nuclear charges and currents. The sum over τ extends from 0 (isoscalar-coupling) to 1 (isovector coupling). Isoscalar and isovector contributions to Eq. (4.4) are proportional to the identity and the third Pauli matrix in isospin space, respectively. For further details, see Appendix B. From Eqs. (2.1) and (2.2), for the -coefficients multiplying S τ , T τ , P τ , Q τ , and R τ we find where s and s label, respectively, the initial and final DM particle spin state, v ⊥ T = v + q/(2µ T ), and µ T is the reduced DM-nucleus mass. Here v is the relative velocity between the incoming DM particle and a reference point in the target nucleus, e.g. its centre of mass. The velocity of individual nucleons with respect to this reference point explicitly appears in the definition of nuclear charges and currents given in Appendix B [12]. With this notation, the squared modulus of T f i admits the following decomposition in terms of nuclear response functions, [12]), and target nucleus spin J [7]:     .2)). Notice also that the interference term proportional to (c τ 10 c τ 18 + c τ 18 c τ 10 ) can only be different from zero when DM couples to quark via scalar and vector mediator simultaneously.

Numerical results
In this section we numerically evaluate the expected rate of DM-nucleus scattering events, Eq. (4.1), assuming xenon as a target material, and focusing on three benchmark models separately: 1) In a first model, h 3 and Im(b 6 ) are the only coupling constants different from zero in Eq. (2.2); 2) In a second model h 3 and Re(b 6 ) are different from zero while all other couplings are zero; 3) In the third model all couplings in Eq. (2.2) are equal to zero, but h 4 and Im(b 7 ). In the three models, DM and the particle mediator have spin 1. When h 3 is different from zero, DM couples to the vector quark current, whereas when h 4 is different from zero DM couples to the axial quark current. Depending on the model, the DM-mediator vertex is determined by the b 6 or b 7 coupling. The first model generates O 19 , the other ones O 20 .
For the three simplified models described above, Fig. 1 Im(b 7 )). . Solid lines correspond to the exact calculation, whereas dashed lines represent the "approximated" results. Interestingly, for both models we find that nuclear recoil energy spectra computed neglecting O 19 or O 20 peak at smaller nuclear recoil energies. Therefore, not only the total number of expected signal events is larger when O 19 and O 20 are properly taken into account, but also the spectral distribution of the events is appreciably modified.

Expected DM relic density
In this section, we briefly review how to compute the DM relic density starting from Eqs. (2.1) and (2.2). We refer to [29,31] (and references therein) for further details on this calculation. Assuming that DM was in kinetic and chemical equilibrium in the early Universe, the time evolution of its cosmological number density, n, is described by the Boltzmann equation,ṅ where H is the Hubble rate in the assumed cosmological model, σ ann is the invariant DM pair annihilation cross section, n eq is the temperature dependent DM equilibrium number density, and v Møl = is the Møller velocity. In Eq. (5.1), angle brackets denote an average over the DM thermal distribution, which can be expressed in terms of the lab frame DM-DM relative velocity, v lab , as follows [31] In the above expression, ε = (s − 4m 2 X )/4m 2 X is the total kinetic energy per unit mass in the lab frame, s = (p 1 + p 2 ) 2 is the squared centre-of-mass energy for the annihilation process, x = m X /T , where T is the CMB temperature, and K 1 (x) and K 2 (x) are the first two modified Bessel functions of the second kind. The product σ ann v lab can be computed as follows where m 3 and m 4 are the masses of the annihilation products, dΩ is an infinitesimal solid angle in the center-of-mass frame, and M is the invariant amplitude for DM pair annihilation. We calculate the invariant amplitude M within the theoretical framework introduced in Sec. 2, starting from Eqs. (2.1) and (2.2). The Boltzmann equation in Eq. (5.1) can be rewritten in terms of the DM abundance Y = n/S, where S is the total entropy density of the Universe, Here, g 1/2 * is the degrees of freedom parameter defined as 6) and the equilibrium density is given by where g is the number of internal degrees of freedom for the DM particle. The quantities h eff and g eff in Eq. (5.6) and (5.7) are effective degrees of freedom for the entropy and energy densities, respectively [31]. An approximate solution to the DM abundance equation, Eq. (5.5), is given by where Y 0 is the present value of the DM abundance (the so-called DM relic abundance), Y f is the DM abundance at the freeze-out temperature, T f , and T f is defined as the temperature below which Y = Y eq . The freeze-out temperature can be found approximately by solving for T [31] 45 π G 1/2 45g where δ ≡ (Y − Y eq )/Y eq is set to 1.5 [31]. Finally, the DM relic density in critical units reads as follows where h = H 0 /100, H 0 is the Hubble constant, and T 0 is the present value of the CMB temperature.

Numerical results
In this section we investigate the compatibility of a signal in a future direct detection experiment, such as DARWIN, with current CMB constraints on the DM relic density in scenarios where the new operators O 19 or O 20 are important. We start by calculating the DM relic density for the simplified model characterised by h 3 and Im(b 6 ) as non zero coupling constants. We perform this calculation setting the model parameters to values producing N = N obs signal events at DARWIN, and keeping the DM particle mass fixed to the benchmark value m X = m * X (which we assume to be reconstructed from the DAR-WIN signal). The constraints N = N obs and m X = m * X (together with the assumption that coupling constants are less than √ 4π, as required by perturbativity) identify a twodimensional surface in the four-dimensional space spanned by the model parameters m X , mediator mass m med = m G , h 3 and Im(b 6 ). If the projection of this two-dimensional surface onto the (m med , Ω DM h 2 ) plane intersects the relic density constraint, Ω DM h 2 = 0.12 [32], then the (h 3 , Im(b 6 )) model is compatible with the observation of N = N obs signal events at DARWIN. Fig. 2 shows the regions in the (m med , Ω DM h 2 ) plane that are compatible with N obs = 1 and m * X = 20 GeV (left panel) and N obs = 150 and m * X = 50 GeV (right panel) when assuming the interval [5 keV, 45 keV] as a signal region and an exposure of 80 ton×year for DARWIN. In both panels, we compare exact results that take into account the O 19 contribution to the observed DARWIN signal (dashed contours) with "approximate" results which neglect the latter (solid contours). We find that neglecting the contribution from O 19 to the expected number of signal events in DARWIN when assessing the compatibility of DARWIN signal and CMB constraints can lead to wrong conclusions, especially for small values of N obs . For example, in the case of N obs = 1 (or, analogously, in a scenario where DARWIN has recorded the first "few" signal events) the range of mediator masses within which DARWIN signal and CMB constraints are compatible varies from a narrow window at small mediator masses (neglecting O 19 ) to a much broader interval (taking O 19 into account).
In this section we presented the results of our numerical calculations focusing on the (h 3 , Im(b 6 )) model. However, we find simular results for model (h 4 , Im(b 7 )), which in the non-relativistic limit generates the O 20 operator. At the same time, in the case of the (h 3 , Re(b 6 )) model, which also generates the O 20 operator in the non-relativistic limit, we find that O 20 can safely be ignored when assessing the compatibility of a DARWIN signal with CMB constraints. At second order in the momentum transfer q and at first order in the transverse relative velocity v ⊥ and in S, we find that four additional non-relativistic operators (compared to the ones already listed in Tab. 1) can be constructed from the Galilean invariants in Eq. (3.1), namely These operators do not appear in the non-relativistic reduction of known simplified models for DM interactions, but they are Galilean invariant and Hermitian, and are therefore not forbidden a priori.

Conclusion
In this work, we focused on the non-relativistic reduction of simplified models for spin 1 DM. We explored two cases separately: in the first case, DM couples to quarks via the exchange of a spin 0 mediator; in the second one, the interactions of DM with quarks are mediated by a spin 1 particle. Our goal was to identify features in the phenomenology of simplified models for DM-quark interactions which are specific to spin 1 DM. We were especially interested in simplified model predictions for DM direct detection and the DM relic density. To calculate the expected rate of nuclear recoils in DM direct detection experiments, for each simplified model we expanded the amplitude for DM-nucleon scattering at second order in the momentum transfer (and at first order in the transverse relative velocity), in analogy with previous works on spin 0 and spin 1/2 DM [11,12]. We then calculated the inverse Fourier transform of the non-relativistic amplitude, summed over all constituent nucleons, and obtained the position-space nuclear potential associated with the underlying DM and mediator model. We finally evaluated the matrix element of the obtained potential between initial and final DM-nucleus states to find the non-relativistic cross section for DM-nucleus scattering. To obtain the DM relic density from a given simplified model, we calculated the cross section for DM pair annihilation into quarks using the formalism in [31].
In the case of spin 1 mediators, we found two DM-nucleon interaction operators arising from the non-relativistic reduction of simplified models for spin 1 DM that were not considered in previous studies and are specific to spin 1 DM. They are quadratic in the momentum transfer and depend on a symmetric combination of the DM particle polarisation vectors. Exploring the phenomenology of the new operators, denoted here by O 19 and O 20 , we found that these can have an important impact on the expected rate of DM-nucleus scattering events at DM direct detection experiments, especially when DM is lighter than 50 GeV, and despite being quadratic in the momentum transfer. This is related to the fact that they arise from simplified models which do not generate momentum transfer independent operators at leading order in the non-relativistic expansion. For example, we found that nuclear recoil energy spectra computed by including O 19 and O 20 or neglecting them differ by up to one order of magnitude for nuclear recoil energies larger than about 20 keV when m X < 50 GeV. Furthermore, we found that the new operators can induce appreciable distortions in the expected nuclear recoil energy spectra, such as a shift in the corresponding peaks.
Focusing on a simplified version of the DARWIN experiment, we investigated under what circumstances a signal at DARWIN can be compatible with constraints from CMB data on the DM relic density in simplified models where O 19 and O 20 are important. We found that, when assessing the compatibility of a signal at DARWIN with CMB constraints, neglecting the contribution from O 19 to the expected number of signal events in DARWIN can lead to inaccurate conclusions, especially when the number of signal events at DARWIN, N obs , is small (e.g. O(1)). For example, for N obs = 1 the range of mediator masses within which DARWIN signal and CMB constraints are compatible varies from a narrow region at small mediator masses, when O 19 is neglected, to a broad region, when O 19 is taken into account.
Finally, we concluded our analysis with a brief discussion on the non-relativistic operators for DM-nucleon interactions which can be built in terms of the symmetric combination of DM polarisation vectors S, but which do not arise from the non-relativistic reduction of simplified models for spin 1 DM. We found that, at second order in the momentum transfer and at first order in the transverse relative velocity and in S, four additional non-relativistic operators can be built from the Galilean invariants in Eq. (3.1). While these additional operators do not appear in the non-relativistic limit of the simplified models introduced in Sec. 2, they are Galilean invariant and Hermitian (see Eq. (6.1)) and therefore not forbidden a priori. In this appendix we explicitly derive the operators O 19 and O 20 from the non-relativistic reduction of Eq. (2.2). The first step in this derivation is to "integrate out" the heavy mediator G µ to obtain an effective Lagrangian that is valid at energies below the mediator mass, m G . This is done by replacing the vector field G µ in Eq. (2.2) with the solution of the equation of motion

Acknowledgments
where G µ = G 0 µ , is the leading order solution in a |q|/m G expansion. Assuming DMmediator and quark-mediator interactions of the general form L int = G µ A µ , where A µ can be any of the DM and quark bilinears appearing in Eq. (2.2), a simple calculation gives the following result Here, we will separately focus on the three A µ bilinear combinations that can contribute to the new operators O 19 and O 20 .

A.1 Mediator interaction via Re(b 6 ) and h 3
In this first scenario, the contact DM-nucleon interactions are described by Eq. (A.2) and the A µ bilinear where the spinor field N describes the nucleon. From Eqs. (A.2) and (A.3), we calculate the amplitude for DM-nucleon scattering, M (1) = 4m X m N M NR , which in the non-relativistic limit is given by Here P = p + p = 2p + q is the sum of incoming and outgoing DM momenta and the vectors ε sµ (e s ) are the four (three) dimensional polarization vectors of the DM field. We also used the non-relativistic reduction of the free nucleon spinor, u N , which at leading order in q gives u r N γ µ u r N 2m N δ r r Kδ r r − 2i(S r r N × q) , (A. 6) where K = k + k = 2k − q is the the sum of the incoming and outgoing momenta of the nucleon. Finally q ν = (E p − E p , p − p), where E p and E p are initial and final DM particle energy, respectively. Expanding the contracted indices in Eq. (A.3), and using Eq. (A.5) and the identity −2m X q ν (ε s * ν ε s µ + ε s * µ ε s ν )δ µ0 = (q · e s )[(P − q) · e s ] + [(P + q) · e s ](q · e s ) ≡ 2q · S s s · P , (A.7) at leading order in q we find Finally, by using the relation v ⊥ = P /(2m X ) − K/(2m N ) we obtain the non-relativistic amplitude Comparing Eq. (A.9) with Tab. 1, we find that the above amplitude can also be written as follows where angle brackets denote expectation values between initial and final DM and nucleon spin and isospin states.
A.2 Mediator interaction via Im(b 6 ) and h 3 In this second scenario, the contact DM-nucleon interactions are described by Eq. (A.2) and the A µ bilinear By replacing the bilinear A µ in Eq. (A.11) in to the effective interaction Lagrangian in Eq. (A.2), for the non-relativistic amplitude for DM-nucleon scattering we find the following expression (A.12) By expanding the contracted indices in the above expression, at leading order in q we obtain the amplitude (q · e s )[(P − q) · e s ]δ r r − 1 2m X (q · e s )[(P + q) · e s ]δ r r − 1 2m N (q · e s )(K · e s )δ r r + 1 2m N (q · e s )(K · e s )δ r r . (A.14) After expressing P in terms of v ⊥ and K, we find that terms proportional to K cancel, and obtain 1 m X q · S s s · q δ r r − 1 m N |q| 2 (S s s X · S r r N ) − (q · S s s X )(q · S r r N ) .
(A.15) Finally, by comparing Eq. (A.9) with Tab. 1, we find that the above amplitude can also be written as The term proportional to O 19 was missing in previous derivations, and can be numerically important.
A.3 Mediator interaction via Im(b 7 ) and h 4 In this last scenario, the contact DM-nucleon interactions are described by Eq. (A.2) and the A µ bilinear The non-relativistic amplitude for DM-nucleon scattering, M NR (defined as above), is thus given by Deriving the above expression, we used the non-relativistic expansion of the axial nucleon spinor bilinear u r N γ µ γ 5 u r N 2S r r N · K 4m N S r r N .
(A. 19) Expanding the contracted indices in Eq. (A.18), we obtain the following non-relativistic amplitude (q · S s s X )(S r r N · K) + 2iq 0 (S s s X · S r r N ) + 1 m X (q × S r r N ) · e s e s · (P + q) − e s (e s · (P − q)) . (A.20) Using the next-to-leading order relation q 0 = P · q/(2m X ) and standard vector identities, Eq. (A.20) reads which can be further simplified by combining second and third term. This leads to the amplitude i m X (q · S s s X )(S r r N · P ) + 2 m X (q × S r r N ) · S s s · q .
(A. 22) Finally, the terms proportional to K and P can be rewritten in terms of a single term proportional to v ⊥ , 2i(q · S s s X )(S r r N · v ⊥ ) + 2 m X (q × S r r N ) · S s s · q . (A.23) Comparison between Eq. (A.23) and Tab. 1 shows that we can rewrite the above Lagrangian as follows Collecting the results found in this section, for the new operators O 19 and O 20 we obtain the following coefficients (A.25)

B From nucleons to nuclei
In this section we provide further details on the relation between differential cross section for DM-nucleus scattering and the underlying nuclear charges and currents (see Eq. (4.4)). We also provide additional information on the DM response functions upon which such cross section depends. The nuclear charges and currents appearing in Eq. (4.4) are defined as follows where x i is the position vector of the i-th nucleon bound in the nucleus, the gradient operator − → ∇ i acts on the coordinates of the x i vector, and σ i is a vector of Pauli matrices acting on the spin space of the i-th constituent nucleon. The above nuclear charges and currents arise from the expansion in spherical harmonics of the nuclear potential that we obtained from M (0) and M (1) , and involve a summation over all constituent nucleons (here