Compatibility of a dark matter discovery at XENONnT/LZ with the WIMP thermal production mechanism

The discovery of dark matter (DM) at XENONnT or LZ would place constraints on DM particle mass and coupling constants. It is interesting to ask when these constraints can be compatible with the DM thermal production mechanism. We address this question within the most general set of renormalisable models that preserve Lorentz and gauge symmetry, and that extend the Standard Model by one DM candidate of mass $m_{\rm DM}$ and one particle of mass $M_{med}$ mediating DM-quark interactions. Our analysis divides into two parts. First, we postulate that XENONnT/LZ has detected $\mu_S\sim\mathcal{O}(100)$ signal events, and use this input to calculate the DM relic density, $\Omega_{DM} h^2$. Then, we identify the regions in the $M_{med} - \Omega_{DM} h^2$ plane which are compatible with the observed signal and with current CMB data. We find that for most of the models considered here, $\mathcal{O}(100)$ signal events at XENONnT/LZ and the DM thermal production are only compatible for resonant DM annihilations, i.e. for $M_{med}\simeq2 m_{DM}$. In this case, XENONnT/LZ would be able to simultaneously measure $m_{DM}$ and $M_{med}$. We also discuss the dependence of our results on $m_{DM}$, $\mu_S$ and the DM spin, and provide analytic expressions for annihilation cross-sections and mediator decay widths for all models considered in this study.


I. INTRODUCTION
The evidence for the presence of dark matter (DM) in the Universe is based upon the observation of anomalous gravitational effects in astronomical and cosmological systems [1]. These systems range from stellar populations in the solar neighborhood [2] to the large scale structure of the Universe [3]. The cosmic microwave background (CMB) radiation is one of the key physical observables in this context [4]. In particular, measurements of the CMB angular power spectrum performed by the Planck satellite set the DM relic density to critical density ratio to Ω DM h 2 ¼ 0.1199 AE 0.0022, with a relative error of less than 2% [5]. A particle in the few GeV up to about 300 TeV mass range with coupling constants at the weak scale is expected to be in thermal equilibrium in the early Universe, and to (chemically) decouple from the thermal bath with a relic density which is typically within a factor of a few from the observed value of Ω DM h 2 [6,7]. This mechanism to generate Ω DM h 2 via chemical decoupling is called DM thermal production. Weakly interacting massive particles (WIMPs) are a typical example of thermally produced DM particles, and by far the most extensively studied candidates for particle DM. For a recent review on WIMPs as DM candidates, see [8,9]. Alternative production mechanisms, e.g., freeze-in and misalignment mechanism, are reviewed in [10,11]. Here, we will focus on thermally produced WIMPs.
The WIMP relic density is obtained through the numerical solution of a collisional Boltzmann equation [12]. The collisional term in this equation takes into account particle physics processes affecting the rate of DM annihilation. These include resonances and co-annihilations [13]. Motivated by increasingly accurate CMB data on Ω DM h 2 [5], predictions for the DM relic density have improved dramatically in recent years and, besides the already mentioned phenomena, now also include, e.g., next-toleading order QCD corrections [14] and the Sommerfeld enhancement of the DM annihilation cross section due to possible DM self-interactions [15]. The impact of a nonstandard cosmological expansion on the DM relic density, e.g. induced by modifications of General Relativity on cosmological scales, has also been studied in detail [16,17]. Computer programmes implementing accurate calculations of the DM relic density include darksusy [18], MicrOMEGAs [19] and MadDM [20]. See, e.g., [21] for a recent comparison of relic density calculations performed using different computer programs.
Among the strategies designed for DM particle identification, the DM direct detection technique will be crucial in the next decade [22]. It searches for nuclear recoil events induced by the scattering of Milky Way DM particles in low-background detectors. See, e.g., [23][24][25] for a review on the current status of direct detection experiments. The discovery of DM at direct detection experiments would constrain the DM particle mass and coupling constants. Importantly, the same parameters are also constrained by CMB data on the DM relic density. Therefore, it is interesting to ask whether the discovery of DM at XENONnT and the DM thermal production mechanism can be compatible. Here we systematically address this question focusing on DM candidates primarily interacting with quarks. Prospects for DM direct detection via DM-lepton interactions have been discussed in, e.g. [26][27][28].
The aim of this work is to determine under what circumstances the DM thermal production mechanism is compatible with the discovery of DM at XENONnT/LZ. In the analysis, we assume that XENONnT/LZ has detected μ S ∼ Oð100Þ signal events. A comparable number of events is expected for DM-nucleus scattering cross sections just below current exclusion limits. We then use this information on μ S as an input to compute the DM relic density. The latter is computed within the most general set of renormalizable models compatible with Lorentz and gauge symmetry extending the standard model by one DM candidate and one particle mediating DM-quark interactions. For most of the models considered here, we find that μ S ∼ Oð100Þ signal events at XENONnT/LZ are compatible with DM thermal production only for a narrow range of mediator masses, i.e. for M med ≃ 2m DM , where M med and m DM are the mediator and DM particle mass, respectively. In this case, a signal at XENONnT/LZ would determine the value of M med and m DM simultaneously, as long as DM is a WIMP. In the analysis, we also relate the compatibility of DM thermal production and DM discovery at XENONnT/LZ to the DM spin and interaction properties. Finally, we discuss the dependence of our results on m DM and μ S . The paper is organized as follows. In Sec. II, we review the set of single-mediator DM models considered in this study. For each model, the relic density is computed as described in Sec. III under the assumption that XENONnT or LZ has detected DM, constraining the model parameters as illustrated in Sec. IV. We present our results in Sec. V, discuss their dependence on m DM and μ S in Sec. VI, and conclude in Sec. VII. For all models considered in this study, analytic expressions for annihilation cross sections and mediator decay widths are listed in the Appendixes.

II. THEORETICAL FRAMEWORK
We calculate DM relic density and DM-nucleus scattering cross sections within the theoretical framework introduced in [29]. This framework consists of "simplified models" extending the standard model by one DM candidate and one particle that mediates the interactions of DM with quarks. This approach is usually sufficient to obtain a good understanding of the DM particle phenomenology, especially in the analysis of DM direct detection experiments and in the calculation of the DM relic density. Simplified models for DM have also extensively been used in the experimental analysis and theoretical interpretation of LHC data, e.g., [30], and to design strategies to extract DM particle spin [31][32][33][34] and particle/antiparticle nature [35,36] from future data. Generically, simplified models are characterized by the nature of the DM particle, which can be a scalar S, a fermion χ or a vector X μ . If DM has spin 1=2, we assume that the field χ is a Dirac rather than Majorana spinor. Independently of its spin, DM is taken to be neutral with respect to the standard model gauge group, but can carry charge under some additional gauge or discrete symmetry group. In this case, the scalar S and the vector X μ must be complex fields. For example, assuming that DM is odd and all other particles are even under a new Z 2 symmetry would guarantee that DM is stable on cosmological time-scales. With the above notation, the Lagrangian of a simplified model for fermionic DM interacting with quarks q through a vector mediator G μ , e.g. a Z 0 -boson, can be written as follows, where a sum over all quark flavors (q ¼ u, d, c, s, b, t) is understood in the terms involving quark bilinears, and λ 3 , λ 4 , h 3 , h 4 are dimensionless constants. Here we assume universal quark-mediator couplings, i.e. the same h 3 and h 4 for all quark flavors. Lagrangians for all models introduced in [29] and considered here are listed in Appendix A. For a given model, we will only consider cases where two couplings are different from zero at the same time.
In the nonrelativistic limit, the models in Appendix A match onto (linear combinations of) nonrelativistic operators for DM-nucleon interactions. These operators are invariant under Hermitian conjugation and Galilean transformations, i.e. constant shifts of particle velocities, but can in certain cases violate T symmetry, as shown in [37]. They have matrix elements between incoming and outgoing DMnucleon states which can be expressed as combinations of the basic invariants under theses symmetries, where v is the DM-nucleon relative velocity, q the threedimensional momentum transfer in the DM-nucleus scattering, and S N and S χ the spin of the DM and of the nucleon, respectively. DM-nucleon interaction operators have systematically been classified in [37], focusing on spin ≤1=2 DM. This classification has subsequently been extended to spin 1 DM in [29]. At linear order in the DM-nucleon relative velocity and at second order in the momentum transfer, 16 independent operators can be constructed, although not all of them appear as leading operators in the nonrelativistic limit of our simplified models for DMquark interactions. The 16 operators are listed in Table I. For each simplified model in Appendix A, the Hamiltonian for DM-nucleon interactions can be expressed as follows where N ¼ p and N ¼ n denotes coupling to protons and neutrons, respectively, the index i ¼ 1; …; 18 (i ≠ 2 and i ≠ 16; see caption in Table I) labels the interaction type, and the coupling constants c ðNÞ i have dimension mass to the power of −2. From Eq. (3), one can calculate cross sections for DM-nucleus scattering as described in Refs. [37,39,40] in detail. Prospects for DM particle detection in this framework have been studied in [31,[41][42][43][44][45][46], and limits on the coupling constants c ðNÞ i have been derived in, e.g., [47][48][49][50][51][52][53].
The coupling constants c ðNÞ i in Eq. (3) are directly related to the coupling constants of the simplified models in Appendix A as illustrated in Table II and described in Ref. [29]. The couplings to nucleons in Table II, h N i , are related to the quark level couplings by nucleon form factors [29,54,55], These values can have large uncertainties. We will briefly discuss how these uncertainties effect our results in Sec. VI.
Only leading contributions to the c ðNÞ i coefficients are considered in Table II. 1 In particular, we neglect momentum-dependent chiral effective field theory corrections [56,57] and renormalization group effects leading to operator evolution [58]. Furthermore, we do not consider charged mediators and models generating the two additional operators introduced in [29],Ô 17 andÔ 18 . The latter arise in the nonrelativistic limit of simplified models with b 6 ¼ ℜðb 6 Þ.

III. DARK MATTER THERMAL PRODUCTION
In the early Universe, the time evolution of the DM space density, n, is described by the Boltzmann equation [12], where H is the expansion rate of the Universe, σ is the invariant cross section for DM annihilation into standard model particles, a dot denotes derivative with respect to the coordinate time, n eq is the DM density at thermal equilibrium, and is the Møller velocity. In Eq. (6), E 1 and E 2 (p 1 and p 2 ) are the energies (four-momenta) of the annihilating DM  [37,38]. Here we adopt the notation introduced in Sec. II. Canonical spinindependent and spin-dependent interactions correspond to the operatorsÔ 1 andÔ 4 , respectively. OperatorÔ 2 is quadratic inv ⊥ andÔ 16 is a linear combination ofÔ 12 andÔ 15 [39]. These are therefore not reported in the table. The operatorsÔ 17 andÔ 18 only arise for spin 1 WIMPs, and S is a symmetric combination of spin 1 WIMP polarization vectors [29]. For simplicity, we omit the nucleon index in the operator definitions.
Some of the coefficients in Table II differ from those in the published version of Ref. [29]. Currently, a revised version of Ref. [29] is in preparation. In the revised version, the new coefficients will agree with Table II and [47].
particles. Angle brackets in Eq. (5) represent an average over the phase-space density F of the DM particles in the initial state, where T is the time-dependent CMB temperature, p 1 and p 2 are three-dimensional momenta and for the phase-space density F we assumed a Boltzmann distribution. Equation (7) is equivalent to the single-integral formula [12] hσv K 1 and K 2 are the first two modified Bessel functions of the second kind, ε ¼ ðs − 4m 2 DM Þ=4m 2 DM is the total kinetic energy per unit mass in the lab frame, and the DM-DM relative speed in the lab frame, v lab , is given by v lab ¼ 2ε 1=2 ð1 þ εÞ 1=2 =ð1 þ 2εÞ. In general, for a two-particle final state, where m 3 and m 4 are the masses of the two particles in the final state, s ¼ ðp 1 þ p 2 Þ 2 , M is the invariant amplitude, dΩ ¼ dϕd cos θ, and θ and ϕ are center-of-mass scattering angle and associated azimuthal angle, respectively. Here we denote the square modulus of M averaged (summed) over initial (final) spin states by jM j 2 . Explicit expressions for jM j 2 for all simplified models in Appendix A are listed in Appendix B.
In terms of DM abundance, Y ¼ n=S , where S is the total entropy density of the Universe, Eq. (5) reads as follows, where G is Newton's gravitational constant, and Y eq is given by In the above expression, g is the number of DM spin states and where g eff and h eff are the effective number of degrees of freedom for the total energy and entropy densities of the  Table I in the proton/neutron basis. In the second column, we omit the index N for simplicity.

Spin 1 DM Coefficient Scalar mediator
Vector mediator Universe, respectively. In the following, we will assume g 1=2 Ã ≃ h eff =g 1=2 eff ≃ 9.5, which is a fairly good approximation for temperatures above the QCD phase transition [12].
Numerical integration of Eq. (5) shows that there is a critical temperature, T f , such that for The critical temperature T f is called freeze-out temperature and can be estimated by solving for x ¼ x f the equation which follows from Eq.
x < x f , and dΔ=dx is a smooth function of x. Here we set the parameter δ to 1.5, a value found comparing results for x f obtained by solving Eq. (12) with estimates based upon Eq. (15) [12]. For x ≥ x f , an approximate solution to Eq. (12) can be found by setting to zero the Y 2 eq term in the right hand side. Indeed, after freeze-out, Y 2 eq is small compared to Y 2 . In terms of T, the approximate solution reads as follows

IV. DARK MATTER DETECTION AT XENONNT/LZ
In this section, we review the equations that characterize DM direct detection at XENONnT/LZ. The aim is to illustrate which DM parameters can be constrained in case of signal detection. These constraints will then determine whether or not DM detection and thermal production are compatible for the simplified models in Appendix A.
The number of observable photoelectrons per DM interaction in the XENONnT/LZ detector is denoted by S1. The expected rate of DM interactions per unit detector mass is [59] where P is a Poisson distribution of mean νðEÞ, and νðEÞ is the number of expected photoelectrons when a nuclear recoil energy E is deposited in the DM-nucleus interaction. The integer n is the number of actually produced photoelectrons in the scattering process. Finally, G is a Gaussian distribution of mean n and variance ffiffiffi n pσ . Here we extract the function νðEÞ from Fig. (13) in [60], set the singlephotoelectron resolution of the XENONnT/LZ photomultipliers,σ, toσ ¼ 0.4, and assume a constant acceptance, ζðS1Þ ≃ 0.4 [60]. The rate of nuclear recoil events per unit detector mass, dR=dE, reads as follows where v min is the minimum DM velocity required to deposit an energy E in the detector. In the experimental analysis, a secondary scintillation signal produced by electrons generated in the DM scattering and drifted to the top of the XENONnT/LZ detector by an electric field is used for background discrimination (i.e. S2 signal). For simplicity, here we neglect the S2 signal-a simplification motivated by the fact that S1 and S2 are anticorrelated [60]. In Eq. (19), v ⊕ is the velocity of the Earth in the galactic rest frame, and ρ DM is the local DM density. Here we set ρ DM ¼ 0.3 GeV cm −3 and assume a Gaussian distribution truncated at the galactic escape velocity v esc ¼ 533 km s −1 for f, the DM velocity distribution in the galactic rest frame boosted to the detector rest frame. These are standard assumptions within the so-called standard halo model, although larger values for the local DM density are favored by astronomical data, e.g. [61][62][63][64][65][66]. For each simplified model in Appendix A, we calculate the differential cross section dσ T ðE; jvjÞ=dE using nuclear response functions implemented in DMFormFactor [39]. We extend the sum in Eq. (19) to the seven most abundant Xenon isotopes, with masses and mass fractions denoted here by m T and ξ T , respectively. Finally, we compute the number of signal events, μ S , integrating Eq. (18) from S1 ¼ 3 to S1 ¼ 70, and multiplying the result by the exposure, ε. For the latter, we assume ε ¼ 20 ton × year. 2 The detection of DM particles at XENONnT/LZ would place constraints on the DM particle mass, m DM , and on the combinations of parameters listed in Table II for the models   2 An analysis extending to S1 ¼ 180 has recently been published by XENON100 [53]. described in Appendix A. Equivalently, it would place constraints on m DM and on the effective mass M eff , where M med denotes one of the mediator masses, m ϕ or m G , and g q and g DM are the coupling constants for the DM-DMmediator andq-q-mediator vertices, respectively. When a simplified model is characterised by the coupling constants g q and g DM and by the leading nonrelativistic operatorÔ i , it will here be denoted byÔ i ðg q ; g DM Þ. For simplicity, from here onwards we will omit the index N in the definition of the nonrelativistic operators. Assuming m DM ¼ 50 GeV and μ S ¼ 100 signal events at XENONnT/LZ, m DM and M eff can be extracted from the data with uncertainties of the order of 20%, e.g. [32]. This is illustrated in Fig. 1 for the case of data simulated from modelÔ 8 ðh 3 ; λ 4 Þ (in this calculation we assume the background model in [32], as described in detail in [60]). The relic density Ω DM h 2 depends on the reconstructed values for m DM and M eff . However, Ω DM h 2 in general depends on M med , g q and g DM separately, and not on their combination M eff . Therefore, an error of about 20% on M eff is negligible compared to the uncertainties on Ω DM h 2 arising from the fact that M med , g q and g DM cannot be constrained independently. Consequently, from here onwards we will assume that a signal at XENONnT/LZ would set m DM and M eff to their true values.

V. COMPATIBILITY OF DIRECT DETECTION AND THERMAL PRODUCTION
In this section, we calculate the DM relic density for the models in Appendix A, showing when it is compatible with the detection of Oð100Þ signal events at XENONnT/LZ. The models in Appendix A are characterized by four parameters: g q , g DM , M med and m DM (as already anticipated, we focus on scenarios where only two coupling constants are different from zero at the same time; see Table II). Here we set m DM and M eff ¼ 0.1M med = ffiffiffiffiffiffiffiffiffiffiffiffi ffi g q g DM p to reference values that we assume to be reconstructed from a hypothetical signal at XENONnT/LZ. Due to parameter degeneracies, it is not possible to associate a unique value of Ω DM h 2 to a ðm DM ; M med Þ pair. Therefore, we compute Ω DM h 2 as follows. We set m DM to its reference value and vary M med in an interval of interest around m DM . We then set g q to minð0.01M 2 med =ðg DM M 2 eff Þ; ffiffiffiffiffi ffi 4π p Þ, as required by the XENONnT/LZ input, and vary g DM in the ð10 −7 ; ffiffiffiffiffi ffi 4π p Þ range, which implies perturbative values for the coupling constants. Furthermore, we impose that the mediator decay width is less than the mediator mass. Through this procedure, we are able to map a signal at XENONnT/ LZ into a region in the M med − x f , M med − hσv Møl i and M med − Ω DM h 2 planes. We perform this calculation for all models in Table II, assuming m DM ¼ 50 GeV as a reference value if not otherwise specified. For M eff , we assume the value that would ideally produce 150 nuclear recoils when g q ¼ g DM ¼ 0.1. We estimate this value by integrating Eq. (19) from 5 to 45 keV in order to allow for a direct comparison of our results with those in [29]. The values of M eff found in this way are listed in Table III for all models  in Table II. These values correspond to μ S ∼ Oð100Þ events at XENONnT/LZ, as one can verify using Eq. (18). In Sec. VI, we will comment on the dependence of our results on the assumed number of signal events. Figure 2 shows x f as a function of M med for different values of g DM and for m DM ¼ 50 GeV. The underlying model isÔ 1 ðh 1 ; g 1 Þ, and g q has been set to the XENONnT/ LZ input. The envelope of the family of curves in the figure identifies a region in the M med − x f plane which is compatible with the observed signal at XENONnT/LZ. All curves peak at mediator masses around 100 GeV, as expected for kinematical reasons: for M med ¼ 100 GeV, the DM annihilation into aqq pair is resonant. We obtain analogous curves in the M med −hσv Møl i and M med − Ω DM h 2 planes. For large M med , the x f envelope tends to a line,  Table III to data simulated from modelÔ 8 ðh 3 ; λ 4 Þ.
In the simulation, we assume m DM ¼ 50 GeV and μ S ¼ 100. When modelÔ 8 ðh 3 ; λ 4 Þ is fitted to the data, the error on M eff is small compared to the uncertainties on Ω DM h 2 arising from the fact that M med , g q and g DM cannot be constrained independently. The figure also shows the bias on the best fit value for M eff arising when a model different fromÔ 8 ðh 3 ; λ 4 Þ is fitted to the data. since in this limit scattering and annihilation cross section depend on the same combination of model parameters (Effective Field Theory limit). Figure 3 compares modelsÔ 1 ðh 1 ; g 1 Þ andÔ 10 ðh 2 ; g 1 Þ in the M med − x f , M med − hσv Møl i and M med − Ω DM h 2 planes. In all panels, we only report the envelopes of the family of curves generated by the g DM parameter. For both models, and in all panels, we observe the expected resonance around M med ¼ 100 GeV, and a plateau corresponding to the Effective Field Theory limit for larger mediator masses. In all panels, we report both exact results and estimates based upon a nonrelativistic approximation for hσv Møl i introduced in Appendix C. For modelÔ 10 ðh 2 ; g 1 Þ, the nonrelativistic approximation is very good, whereas it breaks down in the case of modelÔ 1 ðh 1 ; g 1 Þ due to the presence of a sharp resonance in the annihilation cross section.
The left (right) panel in Fig. 4 shows the x f (Ω DM h 2 ) envelope as a function of M med for the simplified models corresponding to spin 0 DM. In the figure, models are labeled according to the leading nonrelativistic operator for DM-nucleon interactions. ModelsÔ 7 ðh 4 ; g 4 Þ and O 10 ðh 2 ; g 1 Þ are not compatible with the thermal production mechanism for any value of M med , yielding a value for Ω DM h 2 (x f ) much smaller (larger) than the observed one. On the other hand, modelsÔ 1 ðh 1 ; g 1 Þ andÔ 1 ðh 3 ; g 4 Þ generate values for Ω DM h 2 which are in general too large. However, for M med ∼ 100 GeV, i.e. at resonance, direct detection and thermal production can be compatible for modelŝ O 1 ðh 1 ; g 1 Þ andÔ 1 ðh 3 ; g 4 Þ. This result can be interpreted as follows: if DM has spin 0 and is produced thermally, the detection of Oð100Þ signal events at XENONnT would simultaneously determine the DM particle mass, and the mediator mass and spin (at least within the simplified model framework considered in this analysis).
Finally, Fig. 6 shows the regions in the M med − Ω DM h 2 plane corresponding to the detection of Oð100Þ signal events at XENONnT/LZ for simplified models where DM has spin 1. Also in this case, five models can yield TABLE III. Benchmark points ideally producing 150 nuclear recoil events at XENONnT/LZ for m χ ¼ 50 GeV and g q ¼ g DM ¼ 0.1 [32]. We find these values by integrating Eq. (19) from 5 to 45 keV in order to allow for a direct comparison of our results with [29]. The second column shows the leading nonrelativistic operator for the benchmark model. Third and fourth columns report g q and g DM , i.e. the coupling constants for the DM-DM-mediator andq-q-mediator vertices, respectively.

VI. DISCUSSION
In this section, we briefly comment on the dependence of our results on DM particle mass and number of signal events.   uncertainties of up to 30% [54][55][56]. This corresponds to uncertainties in the benchmarks listed in Table III of up to 15%. Varying the input of our numerical simulations for several of the previously discussed models accordingly, we find that these uncertainties will only lead to minor shifts of the predicted regions for DM relic density (less than 1 order of magnitude).

VII. CONCLUSION
We determined under what circumstances the detection of Oð100Þ signal events at XENONnT/LZ can be compatible with the DM thermal production mechanism. The relic density calculation was performed within the most general set of renormalizable models that preserve Lorentz and gauge symmetry, and that extend the standard model by one DM candidate and one particle mediating DM-quark interactions. In this calculation, the detection of Oð100Þ signal events at XENONnT/LZ was used as an input constraining the underlying model parameters.
By increasing the DM particles mass, compatibility regions move towards the bottom right corner in the M med − Ω DM h 2 plane. In the same plane, they move upwards if the number of signal events decreases. Interestingly, for thermal DM models yielding a correct DM relic density only at resonance, direct detection experiments will be able to simultaneously reconstruct DM and mediator mass, in case of signal detection. For these models, M med ≃ 2m DM . Whether or not this value for the mediator mass is excluded by the LHC searches for new physics beyond the standard models is a modeldependent question, which crucially depends on the UV completion of the simplified models discussed in this work. To address this question goes beyond the scope of the present analysis. However, in the large mediator mass limit, M med ≫ m DM , the impact of a XENONnT/LZ signal on future LHC mono-jet searches has recently been studied in [32] for the same set of simplified models considered in this work.
We complemented this analysis by providing analytic expressions for annihilation cross sections and mediator decay widths for all models considered in this study (see Appendix B).

ACKNOWLEDGMENTS
This work has been supported by the Knut and Alice Wallenberg Foundation and is performed in the context of the Swedish Consortium for Direct Detection of Dark Matter (SweDCube). We would like to thank Katherine Freese and Sebastian Baum for interesting discussions on topics related to this work and for their feedback on a first draft of this article. This investigation has also been supported by the Munich Institute for Astroand Particle Physics (MIAPP) within the Deutsche Forschungsgemeinschaft (DFG) cluster of excellence "Origin and Structure of the Universe." Finally, we would like to thank the participants of the program "Astro-, Particle and Nuclear Physics of Dark Matter Direct Detection," hosted by MIAPP in 2017, for many valuable discussions.

APPENDIX A: LAGRANGIANS OF SIMPLIFIED MODELS
In this appendix, we list the Lagrangians considered in our analysis [29].

Scalar dark matter
Scalar and pseudoscalar mediator: Vector and axial-vector mediator: 2. Fermionic dark matter Scalar and pseudoscalar mediator: Vector and axial-vector mediator: 3. Vector dark matter Scalar and pseudoscalar mediator: Vector and axial-vector mediator: APPENDIX B: AMPLITUDES The differential decay width, dΓ, of a particle at rest with mass M into two particles with identical mass m is given by where M is the corresponding amplitude, k is the sum of the incoming four-momenta, and p i and E i , i ¼ 1, 2, are the final state four-momenta and energies, respectively. If M is isotropic, one can integrate Eq. (B1) over d 3 p 1 d 3 p 2 and obtain The relation between DM annihilation cross section in the center-of-mass frame times lab frame relative velocity and the corresponding Feynman amplitude is given in Eq. (10). For all simplified models in Appendix A, we list the modulus squared of the amplitudes for mediator decay and DM annihilation. Expressions can be found in different subsections depending on the DM spin. Results presented in this appendix were cross-checked analytically using FeynCalc [67,68] and through numerical calculations performed with our modified version of WHIZARD [69]. In the amplitudes, we included the averaging over initial state spin and polarization, where applicable. The amplitudes including quarks are for one flavor and without color factors.

Scalar dark matter a. Scalar and pseudoscalar mediator
In the case of scalar DM and scalar or pseudoscalar mediator, the relevant terms in the Lagrangian are −g 1 m S S † Sϕ, −h 1q qϕ and −ih 2q γ 5 qϕ. Assuming that only the coupling constants ðg 1 ; h 1 Þ are different from zero, we obtain and when only the coupling constants ðg 1 ; h 2 Þ are different from zero, we find b. Vector and axial-vector mediator In this case, the relevant terms in the Lagrangian are −ig 4 ðS † ∂ μ S − ∂ μ S † SÞG μ , −h 3 ðqγ μ qÞG μ and −h 4 ðqγ μ γ 5 qÞG μ . For the decay of vector mediators, we have to average over the initial polarization states, which requires the identity When only the coupling constants ðh 3 ; g 4 Þ are different from zero, we obtain When only the coupling constants ðh 4 ; g 4 Þ are different from zero, we find 2. Fermionic dark matter a. Scalar and pseudoscalar mediator In the case of fermionic DM and scalar or pseudoscalar mediator, the relevant terms in the Lagrangian are −λ 1 ϕχχ, −iλ 2 ϕχγ 5 χ, −h 1 ϕqq and −ih 2 ϕqγ 5 q. When only the coupling constants (λ 1 , h 1 ) are different from zero, we find If only the coupling constants (λ 2 , h 2 ) are different from zero, we obtain When only the coupling constants (λ 1 , h 2 ) are different from zero, we find and for λ 2 ≠ 0 and h 1 ≠ 0 we obtain Amplitudes for the process ϕ →qq with h 1 ≠ 0 and h 2 ≠ 0, and for the process ϕ →χχ with λ 1 ≠ 0 and λ 2 ≠ 0 can be found in Eqs. (B4), (B7), (B16) and (B18), respectively.