Dark matter direct detection of a fermionic singlet at one loop

The strong direct detection limits could be pointing to dark matter -- nucleus scattering at loop level. We study in detail the prototype example of an electroweak singlet (Dirac or Majorana) dark matter fermion coupled to an extended dark sector, which is composed of a new fermion and a new scalar. Given the strong limits on colored particles from direct and indirect searches we assume that the fields of the new dark sector are color singlets. We outline the possible simplified models, including the well-motivated cases in which the extra scalar or fermion is a Standard Model particle, as well as the possible connection to neutrino masses. We compute the contributions to direct detection from the photon, the $Z$ and the Higgs penguins for arbitrary quantum numbers of the dark sector. Furthermore, we derive compact expressions in certain limits, i.e., when all new particles are heavier than the dark matter mass and when the fermion running in the loop is light, like a Standard Model lepton. We study in detail the predicted direct detection rate and how current and future direct detection limits constrain the model parameters. In case dark matter couples directly to Standard Model leptons we find an interesting interplay between lepton flavor violation, direct detection and the observed relic abundance.


Introduction
Direct detection (DD) experiments search for dark matter (DM) scatterings off nuclei in underground detectors. The current limits impose very strong constraints on the parameters of Weakly Interacting Massive Particles (WIMPs), which are one of the prototype DM candidates. The current most stringent DD limits for WIMPs in the mass range of [10,1000] GeV come from xenon experiments [1][2][3]. In this work we hypothesize that the absence of DD signals may be reconciled with the WIMP paradigm by generating the scattering at one-loop order and thus with an extra 1/(16π 2 ) 2 suppression of the cross section. As we will see, current and next-generation experiments are able to test significant regions in parameter space of this class of scenarios.
There have been several works in the literature on DD at one-loop order. In Refs. [4][5][6][7] the authors studied DD limits from photon interactions in the context of flavored DM and in Ref. [8] in the context of a radiative neutrino mass model (the scotogenic model [9]) with inelastic Majorana DM. In Ref. [10] the authors performed a detailed study of one-loop scenarios with a charged mediator directly coupled to Standard Model (SM) fields, including the Z and Higgs boson contributions. For couplings to the first and second generation of quarks the dominant contribution may be due to scattering at tree level, while box diagrams may be significant for third generation quarks. Similarly, Ref. [11] studied direct detection of Majorana DM directly coupled to both leftand right-handed SM leptons via two charged scalar mediators. The Z and Higgs contributions were also computed for the scotogenic model in Ref. [12] and also for DM connected to the SM via a neutrino-portal in Ref. [13]. In Ref. [14] the authors studied the one-loop contributions to DD in models with pseudo-scalar mediators or inelastic scattering. In the context of supersymmetry detailed computations have been performed for the bino [15] and wino [16][17][18] DM cases. In the latter scenario loop contributions to DM-nucleus scattering due to gauge bosons may give significant corrections.
In this work we study the DD scattering rate for the case of DM being a SM singlet Dirac or Majorana fermion ψ, which is coupled to a more complex dark sector. A conserved global U (1) or Z 2 symmetry is assumed in order to stabilize the DM particle. In our scenario there are no tree-level contributions to the DD cross section. The lowest order scattering off nuclei occurs at one-loop order via the penguin diagrams in Fig. 1, with a dark fermion F and a dark scalar S running in the loop. We assume that the new particles are color singlets, so that there are no flavor changing neutral currents in the quark sector, and there are only weak limits from direct production at the Large Hadron Collider (LHC). In this way box-diagram contributions to the scattering amplitude are absent. Our main goal is to study analytically the different contributions to the DM-nucleus scattering, as well as to outline possible simplified models, including those with SM fields. In addition we analyze the current limits from DD, as well as constraints coming from the relic abundance, lepton flavor violation (LFV) and anomalous magnetic dipole moments (AMMs). 1 The paper is structured as follows: In Sec. 2 we study the UV completions of the fermionic DM scenario including models with SM particles in the loop. In order to fix the notation we review in Sec. 3 the relevant effective operators for DD at the quark level and also their non-relativistic (NR) versions at the nucleon level. In Sec. 4 we derive analytical expressions for the Wilson coefficients and provide compact expressions in certain limits. In Sec. 5 we perform a numerical analysis of the phenomenology relevant for DD. First we show some numerical examples for the Wilson coefficients at the quark and nucleon level (the latter in their NR version). Afterwards we derive the current limits on the model parameters and discuss future expected sensitivity. We also discuss limits from LFV processes for models in which DM is directly coupled to SM leptons. Sections 4 and 5 contain the main results of this paper. We discuss other phenomenological aspects of the proposed scenario, such as the DM relic abundance, invisible decays and searches at colliders in Sec. 6. Finally we present our conclusions in Sec. 7.
The manuscript also includes several appendices with technical details. The generalization to larger symmetry groups in the dark sector is presented in App. A. In App. B we show a compact expression for the differential cross section in order to make contact with the literature and we Dark scalar Table 1: Particle content and quantum numbers of the fermionic DM scenario with a dark fermion and a dark scalar. The dark sector is charged under a global U(1) symmetry which stabilizes the DM ψ.
briefly review the differential event rate for DD. In App. C we give generic expressions for Higgs and Z boson invisible decays into DM. Relevant formulae for LFV observables and for AMMs of leptons are provided in App. D. Details about the calculation of the relic abundance are collected in App. E and the numerical expressions for the matching to NR operators are given in App. F.

Fermionic singlet dark matter
In the following sections we first present simplified models of Dirac and Majorana fermion DM with vector-like fermions in the loop and then discuss SM particles in the loop.

Dirac dark matter
The new particles can have different combinations of quantum numbers as displayed in Tab. 1. We consider a global U(1) dm symmetry in the dark sector to stabilize DM. It can equally be replaced by a discrete Z n subgroup. Other symmetry groups are discussed in App. A. The interaction Lagrangian for the fields ψ, F and S reads where H is the SM Higgs doublet 2 and V(S, H) denotes the scalar potential. The DM ψ is a SM fermion singlet, but is charged under the dark sector symmetry. The fields F and S are charged under the electroweak gauge group. Electroweak gauge invariance requires them to be in the same SU(2) L irreducible representation of dimension d F and to have equal hypercharge Y F . Notice that in some cases there can be interactions with the SM fields which are subject to strong constraints. We discuss such cases in Sec. (2.3).
In the case of a global symmetry, even if DM is stable at the renormalizable level, higher-order Planck-scale suppressed operators may induce its decay [19]. In particular for a Dirac fermion ψ the dimension-5 operator ψH † ( / DL) with the SM lepton doublet L is one such example. One can construct UV completions of such operators by softly-breaking the global symmetry in the dark sector which induces decays, possibly radiatively. The limits dramatically depend on the DM mass and the Wilson coefficient of the operator. For the rest of the paper we assume that DM is cosmologically stable and that it satisfies all indirect detection constraints on decaying DM.
In our simplified scenario with the interactions given in Eq. (1) it would seem that two of the three new states were stable: ψ and one of S or F . For the following discussion let us assume m S ≥ m F +m ψ so that F is potentially stable, while S can decay. 3 Then, there are two possibilities: (i) If the fermion F is a SM singlet (but charged under the dark group), it also contributes to the DM relic density. 4 Hence, the DD rate of the ψ has to be rescaled by its smaller density under the assumption that the global density scales as the local one, and there is a similar DD rate for F via Higgs penguins with ψ in the loop. (ii) If F is charged under the SM group (SU(2) L charges and/or hypercharge) its electrically charged components have to decay given the stringent limits on charged stable particles [21][22][23][24]. If the components of F mix with SM leptons, they decay like in the model discussed in Refs. [25,26]. Otherwise, as for the DM via the interaction ψH † ( / DL), the fermion F may also decay into SM particles via non-renormalizable operators, which are allowed on general grounds, unless F carries fractional electric charge, or other symmetries forbid them. In this case the fermion F has to decay much faster than the long-lived DM particle. 5 If F is a SM lepton, a charged lepton or a neutrino, it may be stable. Similarly, if F is a right-handed neutrino, it mixes with SM neutrinos and decay. Also, in the case in which S is the SM Higgs doublet and F a heavier fermion, the latter may decay into the SM Higgs boson. We discuss all these possibilities in more detail below.
In general the SM Higgs boson couples to the new scalar multiplet S via a Higgs portal interaction in the scalar potential V(S, H). Depending on the quantum numbers of the particles in the dark sector, it may also have an interaction with the fermion F , for instance if the latter is a SM lepton. In the case of a charged lepton in the loop, the largest Higgs interactions are proportional to the square of its mass (m /v) 2 1, with the electroweak vacuum expectation value (VEV) v 246 GeV. Therefore this contribution is suppressed and can be safely neglected. While these interactions are even further suppressed for Dirac neutrinos, in principle it is possible to have O(1) Yukawa couplings for Majorana neutrinos (we discuss this case in Sec. 2.3.3).
The Higgs portal contribution depends on the coupling of the SM Higgs boson h to a pair of scalars S after electroweak symmetry breaking. In the case of a complex scalar S, we parameterize it in terms of and similarly for a real scalar S with an additional factor 1/2 in order for the Feynman rule (and therefore the expression of the Wilson coefficients) to be identical In the case of a complex scalar S, the Higgs couplings of Eq. (2) are induced by SM gauge invariant Higgs portal interactions such as 3 Similar arguments apply to the other case where S might be stable. 4 In this case, if m ψ mF , coannihilations play an important role [20]. 5 Naively, the scalar S being lighter than the fermion F appears to be more natural given that there are 13 dimension-5 operators which induce decay for a scalar compared to one for a fermion [19]. Table 2: Particle content and quantum numbers of the Majorana DM scenario with a dark fermion and a dark scalar. The dark sector is charged under a Z 2 symmetry which stabilizes the DM ψ.

Dark sector
The term in the first line is always present, while those in the second and third lines require S to be an SU(2) L doublet, S ≡ (S u , S d ) T . Moreover the term in Eq. (6) assumes that S has the same hypercharge as the SM Higgs doublet, which we write after spontaneous electroweak symmetry breaking as H ≡ (0, (h + v)/ √ 2) T . Finally the term in the last line exists for electroweak triplets S ≡ S · σ, where S ± denotes the coefficients of σ ± . 6 In the following we parameterize all the results in terms of λ HS , which allows to easily generalize the result of Higgs penguins for arbitrary combinations of Higgs portals. If S ± (S d,r and S d,i ) have the same mass, their contribution from the interactions (5) and (6) to the DD scattering amplitude exactly cancels due to the relative minus sign in the interaction term. 7 For equal masses the effective coupling λ HS can be generalized from the singlet case to an arbitrary SU(2) L representation of dimension d F by replacing where λ HS,1 (λ HS,2 ) is the coupling of the quartic scalar coupling in Eq. (4) (Eq. (5)).

Majorana dark matter
If the DM particle is in a real representation of a stabilizing dark sector group, it could be a Majorana particle ψ ≡ ψ L + (ψ L ) c (keeping the 4-component notation). We consider the simplest case of a Z 2 symmetry in the dark sector and comment on the general case in App. A. The particle content for Majorana DM is listed in Tab. 2. The Lagrangian is given by If additionally Y F = 0 and consequently S and F both transform according to a real representation, they can be chosen to be a real scalar and a Majorana fermion F = F R + (F R ) c , respectively, and the fermionic part of the Lagrangian simplifies to with y = y 1 = y 2 .
Sector Table 3: Particle content and quantum numbers of the dark fermion scenario with the SM lepton doublet L and a dark scalar.

Standard Model particles in the loop
It is also interesting to study the case where one of the particles in the loop is a SM state. As either the scalar S or the fermion F need to be charged under the dark symmetry, only one of them can be substituted by a SM field. We discuss in the following the cases of S being the Higgs doublet H, and F being the lepton doublet L, the right-handed charged lepton e R or a right-handed neutrino ν R . Interestingly, these types of leptophilic models have some very nice features: (i) the absence of charged stable particles; (ii) the possibility to generate the correct relic abundance by annihilations into leptons; (iii) an interplay with LFV and leptonic AMMs; (iv) the possible relation to lepton number violation (LNV) and neutrino masses; (v) other possible phenomenological signals at future lepton colliders, like MET searches.

Left-handed lepton doublet
The quantum numbers of the remaining states are fixed by demanding that the fermion F in the loop is the SM lepton doublet L, as can be seen in Tab In general, for direct couplings to leptons, it is possible to assign lepton number either to the DM particle ψ or the scalar S. An example with Majorana fermion DM ψ and a discrete Z 2 symmetry (S → −S, ψ → −ψ) is the well-known scotogenic model, proposed in Ref. [9] and extensively studied, e.g., in Refs. [8,[28][29][30][31][32][33][34][35][36][37]. See also the recent review on radiative neutrino mass models [38]. In this case lepton number is broken by the combination of the Majorana mass term of ψ and the operator in Eq. (6). These interactions generate neutrino masses and lepton mixing at one-loop order, which significantly constrain the parameter space of the model. However, in general DD and neutrino masses decouple, because the LNV coupling in the potential could be made arbitrarily small without affecting DD. For fermionic DM, typically, either coannihilations [20] or the freeze-in mechanism [39,40] need to be invoked in order to be compatible with low energy constraints, specially the limit stemming from non-observation of µ → eγ. Table 4: Particle content and quantum numbers of the dark fermion scenario with the SM righthanded charged lepton e R and a dark scalar.

Right-handed charged lepton
If F is the SM right-handed (RH) charged lepton e R , the quantum numbers are fixed as shown in Tab. 4. In this case y 2 = 0 in Eq. (1) for Dirac DM (or eq. (9) for Majorana DM), because the fermions have RH chirality. As in the previous case one should expect new contributions to lepton AMMs and LFV processes. By demanding that the scalar singlet S has lepton number +1, the total lepton number is conserved at the renormalizable level (the term in Eq. (6)) is absent) and consequently no Majorana neutrino masses are induced.

Right-handed neutrino
Dark matter may also couple to right-handed neutrinos ν R with y 2 = 0 in Eq. (1) for Dirac DM (or eq. (9) for Majorana DM). In this case the quantum numbers are fixed as shown in Tab  bosons via the mixing of left-and right-handed neutrinos. This mixing is induced after electroweak symmetry breaking by In this scenario there are two possibilities regarding the nature of neutrinos: they are Dirac fermions for M R = 0, or Majorana fermions for M R = 0. In the latter case, Majorana masses for the active light neutrinos are generated via the seesaw mechanism. In the seesaw scenario the active-sterile mixing angles are tiny, either due to small Yukawa couplings or large right-handed Majorana neutrino masses, and thus the Z penguin contributions and the additional Higgs penguin contributions are extremely small, which agrees with Eq. (19) of Ref. [13]. A possible way-out is to consider an inverse-seesaw scenario, where the suppression needed to have small neutrino masses  originates from a small LNV Majorana mass term, and not from small Yukawa couplings and/or large right-handed Majorana masses.
As DM couples to the SM particles mainly via neutrinos, this is known as the neutrino portal. It has been studied in detail for general heavy SM singlet Dirac and Majorana fermions ν R in Ref. [13] and also in Refs. [41,42].

Higgs doublet
Finally we consider the case of S being the SM Higgs. This fixes the SM quantum numbers of the new particles, which are shown in Tab. 6. This case is qualitatively different, because the neutral component of the electroweak doublet F and the fermion field ψ mix after electroweak symmetry breaking. The lighter of the two neutral mass eigenstates is the DM particle. The Yukawa interactions with the Higgs necessarily induces tree-level contributions to DD via Higgs and Z boson exchange. Although a tree-level contribution exists, DD may still be dominated by the loop-level induced electric or magnetic dipole moments, because they are long-range interactions.

Effective operators for dark matter direct detection
In the following sections we briefly review the effective operators for DM DD. In Sec. 3.1 we show those involving DM interactions with quarks, while in Sec. 3.2 we briefly discuss their NR versions at the nucleon level.

Wilson coefficients at the quark level
Here we review the necessary notation for the effective interactions of the DM with the quarks. The effective Lagrangian at the quark level for a DM fermion ψ is 8 where c q k are the dimensionful Wilson coefficients with the quark q, c g andc g are the Wilson coefficients for gluon operators and µ ψ and d ψ magnetic and electric dipole moments. We implicitly assume that the operators are generated at a scale above the nuclear scale, ∼ 2 GeV. See App. F for further details.
We focus on the contributions to spin-independent (SI) and spin-dependent (SD) operators of photon, Z boson and Higgs penguins which are not momentum or velocity suppressed. The latter would yield very small rates, as there is already the one-loop squared factor at cross section level, 1/(16π 2 ) 2 . We start the discussion with the case of ψ being a Dirac fermion and later on discuss the case of DM being a Majorana particle.
For SI scattering the relevant dimension-6 effective operators are where q denotes the quark field. O q SS is generated by the gauge-invariant dimension-7 operators (ψψ)(Q LH u R ) and (ψψ)(Q L Hd R ), where Q L , u R , d R represent the quark flavor eigenstates. O q SS flips chirality and it is generated by Higgs exchange and thus we factor out the quark mass m q . O q VV preserves chirality and is generated by photon or Z exchange. The contribution from the photon penguin can be related to the anapole moment ψγ µ ψ ∂ ν F µν and the (non-gauge invariant) milli-charge operator ψγ µ ψ A µ via the equations of motion for the photon.
There are also scatterings of the DM with gluons at two-loop order which generate the dimension-7 operators: where a = 1, ..., 8 are the adjoint color indices, α s is the strong coupling constant, G µν the gluon field strength tensor andG µν ≡ 1 2 µνρσ G ρσ its dual. O g is induced from O q SS after integrating out the heavy quarks. We explicitly factorized out a loop factor, as these operators can never be generated at tree level.
For SD interactions the relevant dimension-6 effective operators are where Only the Z boson contributes to O q AA . In SM effective theory the tensor operator may arise from one of the dimension-7 operators (ψσ µν ψ)(Q LH σ µν u R ) and (ψσ µν ψ)(Q L Hσ µν d R ) which are however not induced at leading order.
Photon penguins also generate long-range interactions which are described by the magnetic (CP-even) and electric (CP-odd) dipole moments of the DM ψ, namely with µ ψ and d ψ the coefficients of the magnetic and electric dipole moment operators introduced in Eq. (12), respectively. The latter are generated radiatively and therefore it is convenient to factorize a loop factor. In the case of a Majorana DM particle there are only operators with the bilinears ψψ, ψ γ 5 ψ and ψ γ µ γ 5 ψ, so that the vector O q VV , the tensor O q TT and the dipole moment operators, O mag and O edm , vanish identically. Thus, for SI scattering only the Higgs penguin which generates O q SS is present. For SD scattering O q AA generated by the Z boson can also be non-vanishing. In this case we also compute the photonic contribution to the anapole operator which gives rise to momentum-suppressed and velocity-suppressed NR operators (both SI and SD). See also Ref. [44] for a study of the phenomenology of Majorana DM in EFT.
In general the penguin contributions are isospin-violating, i.e., with different couplings to protons and neutrons (f n = f p ). This isospin violation is maximal for photon contributions which only couple to protons. The latter dominate the DM-nucleus scattering via the dipole moments µ ψ and d ψ . Hence for SI DM-nucleus scattering the enhancement due to coherent scattering is Z 2 instead of A 2 , with Z (A) being the number of protons (nucleons) of the nucleus.

Non-relativistic Wilson coefficients at the nucleon level
The previous Wilson coefficients at the quark level generate non-trivial Wilson coefficients at the nucleon level [45][46][47]. The different contributions generally interfere. The matrix elements of DMnucleon scattering can be written as a linear combination of the following relevant NR operators in the convention of Ref. [48]. I ψ (I N ) denotes the identity operators for DM (nucleons), S ψ ( S N ) denotes DM (nucleon) spin, and q and v ⊥ describe the momentum and velocity exchange. We use DirectDM [48] to match the simplified models onto the NR operators. The numerical expressions for the matching to NR operators are collected in App. F. The NR Wilson coefficients may depend on the transferred momentum q. Note the different normalizations of the spinors and the effective operators between Refs. [45][46][47]49] and Refs. [48,[50][51][52]. In addition to the different definitions of the quark-and nucleon-level operators, in order to translate between these conventions one needs to multiply the NR Wilson coefficients of Refs. [48,[50][51][52] by 4 m ψ m N (4 m ψ | q| 2 ) in the case of contact (long-range) interactions. Further details can be found in the recent Refs. [48][49][50][51][52]. The differential cross section for DM scattering off nuclei is given in App. B.

Analytical results
The effective operators in Eq. (12) are generated at one-loop order from penguin diagrams mediated by the photon and the Z and Higgs bosons. We have computed the different contributions using the Mathematica packages FeynRules [53], FeynArts [54], FormCalc and LoopTools [55][56][57], ANT [58] and Package X [59,60]. As we show below, although the long-range interactions are expected to dominate, the short-range effective operators become relevant in some cases. One obvious example is DM-nucleus scattering of Majorana DM, since the dipole moments vanish. Therefore we show below all relevant contributions. The interesting SI (SD) interactions in Eq. (12) are given by the dipole moment operators O mag and O edm as well as the operators O q SS , O g and O q VV (O q AA ). All the other operators in Eq. (12) are suppressed in the limit of small momentum transfer by a factor | q| 2 /m 2 N or | q| 2 /m 2 ψ , where m N is the nucleon mass. In the following we express the SI and SD Wilson coefficients in Eq. (12) in terms of the ratios and the loop function It is convenient to define the vector and axial Yukawa couplings: Similarly, the interaction of the fermion F with the Z boson in Eq. (1) may be written in terms of vector and axial-vector couplings, namely where e > 0 is the proton electric charge, and s w (c w ) denotes the sine (cosine) of the weak mixing angle. If F is a vector-like fermion, then we have where Q is the electric charge of the (component of the) field F , in units of e, and Y F is the corresponding hypercharge. Conversely, for a SM lepton F we have and the Yukawa couplings are For simplicity of notation we report the full analytic results for SU(2) L singlets F and S. In the case of no mass splittings between the components of the SU(2) L multiplets of dimension d F it is straightforward to generalize the results: The expressions for photon penguins and electric and magnetic dipole moments are generalized by replacing Q → d F Y F . Higgs penguins are generalized for different scalar multiplets as in Eq. (8).
Most Z penguin contributions (apart from some with chiral SM fermions) vanish at leading order. This is also the case for other SU(2) L multiplets.
We summarize below the relevant contributions to the (Dirac or Majorana) DM-quark scattering amplitude. We have checked that our expressions agree with those reported in the literature in the appropriate limits: dipole and anapole moments in Refs. [5,6,8,10], and also for the Z boson contributions in Refs. [10,15].

Dirac dark matter
The leading contributions for Dirac fermion DM are from dipole moments, the operators O q V V and O q AA , and the scalar operator O q SS . Integrating out heavy quarks induces the gluon operator O g .

Electromagnetic dipole moments
The magnetic and electric dipole moments are given by and Both O mag and O edm flip chirality and therefore the dominant contributions to their coefficients are proportional to the heaviest fermion mass, either m ψ or m F . In the limit m ψ m F < m S these expressions reduce to

Photon penguin
Photon penguins induce the operator O q V V . The relevant Wilson coefficient in the effective Lagrangian (12) is where Q q is the electric charge of the quark q in units of e > 0. In the limit m ψ m F < m S the expression above reduces to In the case the mass of the fermion in the loop is much smaller than the momentum transfer,

Z penguin
For a vector-like fermion the resulting SI and SD scattering amplitudes are suppressed by | q| 2 /m 2 F and | q| 2 /m 2 S due to a cancellation between the diagrams where the Z boson couples to the scalar and to the fermion. Therefore, no strong constraints on the model parameters can be obtained. For SM leptons in the loop we distinguish two cases: (i) If S is a singlet under SU(2) L , the axial-vector coupling in Eq. (27) is z A = 0 and both SI and SD scattering amplitudes are suppressed as for a vector-like fermion.
(ii) If S ≡ (S 0 , S − ) T is a doublet under SU(2) L , there are contributions from both diagrams where the Z boson is attached to the SM lepton or the scalar in the loop The couplings q V,A are q V /e = 3 − 8s 2 w /(12c w s w ) and q A /e = −1/(4c w s w ) for up-type quarks and q V /e = −3 + 4s 2 w /(12c w s w ) and q A /e = 1/(4c w s w ) for down-type quarks. Q f denotes the electric charge of the lepton. We define x ψ ≡ m ψ /m S − , x ≡ m /m S − and x ν ≡ m ν /m S 0 with the charged lepton mass m and the neutrino mass m ν . This agrees with the expression in Ref. [10]. The contribution with light active neutrinos in the loop is negligible because it is proportional to x 2 ν and thus the contribution is entirely determined by the charged lepton in the loop. However, for models with a neutrino portal as outlined in Sec. 2.3.3 there may be a sizable contribution from right-handed neutrinos (mixed with left-handed neutrinos) in the loop. In the limit of small DM mass, x ψ 1, the contribution of right-handed neutrinos is with the active-sterile mixing angle θ. We define x N ≡ m N /m S with the heavy neutrino mass m N . These interactions are also discussed in Ref. [13] (see also Ref. [61]). In the case of mixing of vector-like charged fermions with SM charged leptons, there is an overall minus sign in the expressions of Eqs. (38) and (39).

Higgs penguin
At leading order in | q| 2 there is only the contribution to the SI scattering amplitude. The relevant Wilson coefficient generated by the Higgs-portal interaction is As previously mentioned we neglect the contribution from the Higgs penguin where the Higgs boson couples to a SM lepton in the loop, because it is suppressed by (m /v) 2 1. 9 As in the case of O mag and O edm , the operator O q SS flips chirality, and therefore the dominant contribution to its Wilson coefficients is proportional to either m ψ or m F . If both F and ψ are charged under U(1) dm , then m ψ < m F and thus the largest contribution comes with the chirality flip on the fermion line of F . On the contrary if F is a SM lepton the largest Higgs contribution is proportional to m ψ . In the limit m ψ m F < m S , Eq. (40) simplifies to

Photon penguin
For Majorana DM ψ the electromagnetic dipole moments identically vanish and the only allowed electromagnetic form factor is the anapole moment. This gives rise to the effective operator O q AV in Eq. (17). We obtain In the limit m ψ m F < m S this simplifies to In the case the mass of the fermion in the loop is much smaller than the momentum transfer, m F −q 2 , we have where x q ≡ −q 2 /m 2 S .

Z and Higgs penguin
For a vector-like fermion F in the loop the Z penguin diagram does not contribute to the SI scattering amplitude, because the DM vector current identically vanishes for Majorana fermions. The SD scattering amplitude is suppressed by | q| 2 /m 2 F and | q| 2 /m 2 S due to a cancellation similar to that occurring in the case of Dirac fermion DM, see Sec. 4.1.3. If F is a left-handed lepton doublet, and consequently S ≡ (S 0 , S − ) T is an SU(2) L doublet, we find at leading order in | q| 2 : c q VV = 0 and c q AA is a factor of two larger than result for the Dirac case provided in Eq. (37). If F is a right-handed charged lepton or a right-handed neutrino, the scalar S is necessarily an SU(2) Lsinglet and thus the axial-vector coupling z A in Eq. (27) is zero and both SI and SD scattering Z-mediated amplitudes are suppressed as for a vector-like fermion. In some models, like with right-handed neutrinos or with vector-like fermions, there can be mixing with SM leptons. These generate couplings to the Z and the Higgs bosons, see discussion around Eqs. (38) and (39), and footnote 9.
For the Higgs penguin there is only a contribution to the SI amplitude c q SS at leading order in | q| 2 , which again is a factor of two larger than in the Dirac DM case, given in Eq. (40). The fact that the h and the Z penguin contributions to the non-zero Wilson coefficients, c q SS and c q AA , are a factor of 2 larger for Majorana than for Dirac DM, can be understood from the presence of extra crossed diagrams for Majorana particles, where the initial and final DM particles are interchanged.

Numerical analysis
We use LikeDM [62,63] to compute the differential rates and the experimental upper bounds on our scenarios. We have also performed cross checks with the program of Ref. [49]. First we show results for the event rates and upper limits for Dirac and Majorana DM, having either vector-like fermions or SM leptons in the loop. For the latter case we also show upper limits from LFV signals. In the following we parameterize the vector and axial Yukawa couplings of Eq. (24) in terms of their absolute value and phase as y V = |y V |e iφ V and y A = |y A |e iφ A .

Wilson coefficients at the quark level
In order to illustrate the relative weight of the different contributions, we plot in Fig. 2 10 for which we use E R = 8.59 keV (which is a reasonable value for xenon nuclei, with mass m Xe 132 GeV). The magnetic dipole moment dominates, followed by the photon short-range contribution which is roughly ∼ 10 −8 fm. The increase of µ ψ and the Higgs contribution with m ψ is easily understood from chirality arguments. This also implies that µ ψ and c u SS are suppressed with respect to the case of vector-like leptons (cf. upper-left panel of Fig. 2) by the DM mass, except in the region of m ψ close to m S . The Higgs and the Z penguin interactions are always very suppressed (for the Z penguin the SD amplitude is smaller than the SI contribution, due to the factors q V,A /e in Eqs. (36) and (37)), below 10 −11 fm, and therefore they can be safely neglected. All Wilson coefficients increase for m ψ ∼ m F + m S .
For Majorana DM with SM lepton doublets in the loop (bottom right) the Higgs and the Z SD amplitudes are a factor of two larger than in the Dirac case and with the same dependence on m ψ , while the anapole Wilson coefficient (dot-dashed purple) is slightly larger than the photon short-range contribution present in the Dirac case. Notice that this is the opposite behavior of the case with vector-like fermions.

Wilson coefficients at the nucleon level
The previous Wilson coefficients at the quark level can interfere and generate non-trivial effective operators at the nucleon level, see Sec. 3.2. We plot in Fig. 3          they are generated dominantly by µ ψ , i.e., chirality needs to be violated. Similarly c N 4 (blue) with p (dotted) and n (dashed) have a dominant contribution from µ ψ and therefore increase with m ψ . Regarding c n 1 (dashed black), the increase of its slope reflects the fact that the short-range Higgs contribution (which grows with m ψ ) increasingly becomes more and more comparable to the photon short-range coefficient, but in any case c n 1 remains very suppressed. c p 1 is dominated by the photon penguin, and both the short-range contribution parameterized by c q V V,γ and the longrange contribution from the magnetic moment µ ψ /m ψ are important. Due the dependence of the quark-level Wilson coefficients on the DM mass, the NR Wilson coefficient c p 1 is basically constant with respect to it. For Majorana DM with SM leptons (bottom right) the Wilson coefficients c p 8 (solid purple), c p 9 (solid magenta) and c n 9 (dashed magenta), which are generated by the anapole operator, dominate. c N 6 (c p 6 and c n 6 are superimposed in the plot) do not increase with the DM mass, unlike in the Dirac case, because here they come from c N AA and not from µ ψ ; c N 1 , generated by the Higgs penguin diagram, increases with m ψ and is very similar for n and p (c n 1 and c p 1 are superimposed in the plot). Finally, c N 4 , generated by the Z penguin, are similar for both n and p    (superimposed in the plot) and very suppressed, as expected.

Direct detection event rates
The different Wilson coefficients are expected to generate different features in the DD differential spectrum. In the upper-left panel of Fig. 4 we plot the DD differential event rates in xenon versus the recoil energy E R for Dirac DM with a vector-like fermion F (solid blue) and with a righthanded tau (dotted green), and for Majorana DM with a vector-like fermion (dashed red) and with a right-handed tau lepton (dot-dashed purple). For details on the astrophysical assumptions used in the numerical analysis see Ref. [63]. The rate for Dirac DM with a vector-like fermion is roughly 9 orders of magnitude larger than that with a SM lepton (a tau lepton in this case), because in the latter case there is no contribution from the electric dipole moment d ψ and the magnetic dipole moment µ ψ is suppressed by the DM mass m ψ . The smallest rate occurs for Majorana DM with a right-handed tau in the loop. The relative size of the spectra is obvious from   In the bottom panel of Fig. 4 we plot the DD differential rates for Dirac DM with coupling to a right-handed electron (solid blue), muon (dashed red) and tau (dotted green). The spectra are the largest for the electron (the lightest lepton), with maxima at roughly the same recoil energy. The maxima go approximately in the ratios ∼ (4 : 2 : 1) for e, µ, τ . This is due to the dependence on the short-range contribution of the photon penguin, c q V V,γ , via the NR Wilson coefficient c p 1 . The spectra are dominated by the photon short-range contribution c q V V,γ for this choice of parameters.

Direct detection limits
Next we study the upper limits that current DD experiments can impose on the scenarios discussed so far. In order to illustrate current direct detection limits, we consider different scenarios of TeVscale dark sectors. We also discuss how the limits vary with the masses of the particles in the loop. We show the 90% C.L. upper limits from current DD experiments that have xenon as a target, which provide the most stringent limits for SI interactions for our range of DM masses. We show m ψ 5 GeV, as very light DM does not produce recoils at energies above the threshold of the DD experiments. The limits are subject to large uncertainties from nuclear physics and astrophysics as well as from experimental uncertainties. In the following we do not show limits from Higgs and Z boson invisible decay widths into DM, as those are weaker than the ones coming from DD in our scenarios. In Sec. 6.2 we discuss some examples where these limits can be relevant, and complementary to DD, specially for light DM masses, and in App. C we provide the relevant expressions for the Higgs and Z boson invisible decay rates.
In the left panel of Fig. 5 we plot the upper limits for Dirac DM in the plane |y V | versus m ψ for XENON1T (solid brown), PandaX (dashed green) and LUX (dotted purple), together with their combined limit (thicker solid red line). We have fixed m S = 500 GeV, m F = 600 GeV, λ HS = 0.1, the ratio of Yukawa couplings |y A |/|y V | = 1.3 and the phases of the Yukawa couplings φ V = 0 and φ A = 1.4. As expected the bounds are weakened at very large and very small DM masses. At large DM masses the limits appear to approach a constant value, instead of decreasing as 1/m ψ as expected from the DM number density. This is due to the non-trivial dependence of the Wilson coefficients on m ψ . In particular the Wilson coefficients generally increase for m ψ → m F + m S . The |y V | limits are of the order of ∼ 10 −2 for a large range of DM masses between 10 GeV and 500 GeV. This is a clear example of the superb sensitivity achieved by DD experiments, which are able to probe such small Yukawa couplings for loop-induced scenarios of Dirac DM.
In the right panel of Fig. 5 we show the limits for Majorana DM with λ HS = 0.1 (dashed red) and λ HS = 3 (dotted green), together with those for Dirac DM (solid blue). The current limits for Majorana DM are very weak, close to the naive perturbativity limit. Notice that the Higgs interactions are non-negligible: changing λ HS = 0.1 to λ HS = 3 the upper bound on the Yukawa couplings improves by a factor of ∼ 6 (at the level of the rate, the scalar quartic coupling enters quadratically, while the Yukawa couplings enter to the fourth power). The difference with respect to the strong limits for Dirac DM stems, of course, from the absence of dipole moments for Majorana DM. In the gray shaded region, the Yukawa coupling is non-perturbative, |y V | > √ 4π, and therefore the one-loop computation cannot be trusted.

Interplay with lepton flavor violation and relic abundance
When there are SM charged leptons running in the loop, there may also be limits from LFV processes. We provide the relevant expressions for α → β γ, µ − e conversion and α → β γ δ in App. D. 11 It is therefore interesting to study the interplay between both types of signals. Although one may naively expect that LFV limits are stronger (because an accidental symmetry of the SM is violated), we see in the following that this is not the case in all scenarios.
In Fig. 6, top-left panel, we plot the DD upper limits 12 in the plane |y 1 | versus m ψ , assuming equal couplings to all leptons, i.e., y e 1 = y µ 1 = y τ 1 = y 1 (we denote this the "symmetric" case). In Fig. 6 top-right, middle-left and middle-right panels we show the cases of no couplings to taus, electrons and muons, respectively. Left-handed and right-handed leptons in the loop lead to the same result. We show the cases of Dirac DM (solid red), and Majorana DM with λ HS = 0.1 (dashed light blue) and λ HS = 3 (dashed green). We have fixed m S = 1000 GeV for the four upper plots. The most relevant 90% C.L. LFV limits are shown using dotted lines: µ → eγ (green), µ − e conversion (orange), µ → 3e (black) and ∆a µ (brown). 13 Notice that LFV limits do not depend on whether the DM is a Dirac or Majorana fermion. Also, we emphasize once more that DD limits are subject to large astrophysical and nuclear uncertainties, which are absent in       the case of LFV experiments. In addition we plot the contour of the DM relic abundance, set by t-channel DM annihilations ψψ → α β mediated by the scalar S, with a dot-dashed navy blue (purple) line for Dirac (Majorana) DM, whose leading contribution is from s-wave (p-wave) scattering. We use the instantaneous freeze-out approximation which is sufficient for our purposes (see Sec. 6.3.1 and App. E for more details and the relevant expressions). Above the Ωh 2 contour the DM would be under abundant and requires an additional component of DM to account for the observed relic abundance. Below the Ωh 2 contour DM is over abundant if its abundance is solely set by freeze-out, and thus there has to be a mechanism to further deplete its density. It could be reduced via co-annihilation and resonant effects [20], multi-body scatterings [64][65][66][67][68], or a non-trivial thermal evolution in the early universe [69]. In case ψ does not account for all of the DM abundance the DD limits have to be rescaled appropriately. Assuming thermal freeze-out reproducing the correct relic abundance imposes a lower bound on the DM mass. In the case of equal couplings to all leptons with m S = 1000 GeV, m ψ The main changes in the case of no couplings to taus, electrons and muons (top-right, middleleft and middle-right panels in Fig. 6) are in the LFV limits, as depending on the flavor structure, different processes are possible. In these panels the relic abundance contours are almost identical, as the SM leptons are always much lighter than the DM (and therefore phase space plays no significant role). Of course, the contours are at somewhat larger Yukawa couplings than for the "symmetric" scenario, as in the latter there were more available annihilation channels. The DD limits are also slightly modified due to the different masses of SM leptons in the loop (see also the lower panel of Fig. 4). When there are no couplings to taus (top-right panel), the LFV limits are almost identical to the "symmetric" scenario, because they are driven by the first family. However, for no couplings to electrons or muons (middle panels), DD limits are more stringent than LFV limits for Dirac DM with a mass above a certain value. This is quite remarkable: DD experiments are able to better constrain scenarios where an accidental symmetry of the SM is violated than experiments directly searching for it. Interestingly, limits on |y 1 | from trilepton τ decays (τ → 3 ) dominate over radiative τ decays (τ → γ) in contrast to the limits from muon decays. As the limits from τ decays are generally weaker and thus the corresponding Yukawa couplings larger, box-diagram contributions to trilepton decays may give a sizable contribution and thus break the dipole dominance.
A few interesting remarks can be drawn from these plots. First, note that the DD limits with SM leptons in the loop, even for Dirac DM, are much weaker than in the scenario with vector-like fermions in the loop, as also demonstrated in the top-left panel of Fig. 4. Second, clearly the LFV limits are the strongest ones, with µ → eγ the most stringent among them. Its limit on the Yukawa coupling |y 1 | is a factor of a few stronger than the one of DD for Dirac DM. Again, the DD limits become very strong close to m ψ → m S + m F as in the case with vector-like fermions. Third, for scalar masses at the TeV scale, the DD limit already excludes the production via thermal freeze-out for Majorana DM, and also for Dirac DM in the mass range 5 GeV m ψ 200 GeV. Finally, the muon AMM constraint is always very weak, being the limit above the perturbativity bound.
In Fig. 6, bottom panels, we show two examples of a scalar S in the loop with a different mass: m S = 300 GeV (left plot) and m S = 5000 GeV (right plot). All limits are generically stronger for m S = 300 GeV and weaker for m S = 5000 GeV compared to m S = 1000 GeV. In particular, the relative contribution of the box diagrams and the dipole moment for the trilepton τ decay changes: for m S = 300 GeV τ → eγ sets a stronger limit than τ → 3e. Similarly the Yukawa coupling required to explain the observed relic abundance also has to be larger for heavier scalar masses, as already discussed. Indeed, for m S = 5000 GeV almost all the limits on the Yukawa couplings are in the non-perturbative region.
In summary, strong limits can be set for Dirac DM with vector-like fermions in the loop. For Dirac DM with SM leptons in the loop LFV limits or DD limits may set the strongest bounds depending on the flavor structure and the DM mass. Therefore, the two limits are complementary: LFV limits are more important for DM coupling to both muons and electrons, whereas DD limits dominate if there are no LFV processes of type µ → eX, X being anything, and the DM mass is not too small (m ψ 5 GeV). For Majorana DM, LFV limits, if present, are generally more stringent than constraints from DD. Future DD experiments and LFV limits on τ decays are expected to improve by 1-2 orders of magnitude and hence the situation is not expected to change dramatically. If µ − e conversion in nuclei and/or µ → 3e expected sensitivities (by several orders of magnitude) are achieved, LFV limits will continue to dominate and even increase their difference with respect to DD.

.1 LHC searches
Generally colliders may only set competitive limits via missing energy searches for light DM and SD interactions. In the scenarios discussed here, naively the production of DM particles at the LHC occurs at one-loop level via the penguin diagrams in Fig. 1 and is therefore suppressed. For example, Ref. [25] showed that there are only very weak collider limits on a model with a magnetic moment interaction. Thus it is more promising to search for the mediators S and F at colliders via qq → FF , SS * mediated by the photon, the Z boson and/or the Higgs. If the new fermion and scalar have electric charge, the production is dominated by the Drell-Yan process. Higgs-mediated production of exotic particles has been discussed in e.g. Ref. [70]. As we are assuming that the new particles are not colored, only modest lower limits (below 1 TeV) are expected, unless very large SM quantum numbers (for instance electric charges) are invoked. The dark sector particles may decay invisibly into DM and a lighter dark sector state. The phenomenology of these decays are however model-dependent, see discussion in Sec. 2. Another interesting option would be to search for DM in models with electrons/muons running in the loop at future lepton colliders. The main production process is via t-channel exchange of the scalar, + − → ψψ with = e, µ.

Z and Higgs boson invisible decays
If the DM ψ is sufficiently light [m ψ < m H /2 (m Z /2)] there is an additional contribution to the invisible width of the Higgs (Z) boson. In App. C we present the relevant expressions for these processes. We find that there are no limits from Z or Higgs boson decays into DM for the   [25,53] GeV ( [22,55] GeV) by ATLAS (CMS) for Dirac DM and [16,57] GeV ( [14,58] GeV) by ATLAS (CMS) for Majorana DM.
To summarize, while for vector-like fermions there are no limits, for SM particles in the loop there may be interesting constraints in the absence of LFV. Indeed, there is a well-known complementarity between invisible decays and DD. The experimental energy threshold of DD experiments limits their ability to impose limits for arbitrarily low DM masses and thus invisible decays may set competitive limits for low DM masses.

Relic abundance
The production of the correct relic DM density in the early universe is generally model-dependent. Although it is not the main focus of this work, we briefly outline different avenues to obtain the correct relic density. See e.g. Ref. [74] for a connection of DD with thermal freeze-out.

Thermal freeze-out
If m ψ > m S , m F (but of course m ψ < m S + m F ), the relic abundance can be set via the t-channel interactions ψψ → SS * or ψψ → FF . Subsequently, S and F can decay to SM particles, in some cases at loop level or via non-renormalizable operators. In particular if F is a SM lepton α , DM annihilations to SM leptons ψψ → α¯ β may set the relic abundance. For Dirac DM the cross section is not velocity suppressed and thus the leading (s-wave) part of the thermally averaged annihilation cross section 14 is given by where we have summed over all possible final state leptons (neutrinos and charged leptons) in the limit of vanishing lepton masses. Here i = 1 (2) for couplings to LH (RH) leptons, see Eq. (1). For Majorana DM the annihilation cross section is velocity suppressed and the leading contribution is due to p-wave scattering 15 where x = m ψ /T . As discussed in App. E, for DM masses in the range 10 GeV m ψ 10 4 GeV we obtain the correct relic abundance for cross sections σv | D 14 The thermally averaged cross section σv = a + 6b/x with x = m ψ /T is obtained by integrating over the annihilation cross section σv = a + bv 2 , after it has been expanded up to second order in the relative of velocity of the two DM particles in the center of mass frame v = | v|. Note that, although the DM is non-relativistic at freeze-out, the relative velocity is not small, v f = 12/x f 0.7 c in terms of the speed of light c. 15 Annihilation channels with 3-body final states which lift the velocity suppression are generally not important during freeze-out due to the additional phase space suppression, but they are very important for indirect detection. Their importance for indirect detection has been pointed out in several papers [75,76], see also Refs. [6,77].
If ψ is the lightest particle in the dark sector (i.e., m ψ < m S , m F ), DM may annihilate at oneloop order into quarks via the penguin diagrams in Fig. 1. However this is very suppressed and results in an over abundance of DM and requires another mechanism: (i) In a larger dark sector DM may annihilate into other lighter dark particles, ψψ → XX which subsequently decay to SM particles. These new light particles may lead to large DM self-interactions, see for instance Ref. [78].
(ii) Co-annihilation and resonant effects [20] may increase the effective thermal annihilation cross section. For example processes like ψF → S * → HH with (m F − m ψ )/m ψ 1/20 could be induced by a coupling of S to the SM Higgs. 16 Similarly there may be coannihilations with S. If S has gauge interactions the dominant channel may be SS → SM SM (see for instance Ref. [79]) if (m F − m ψ )/m ψ 1/20 and ψ and S are in thermal equilibrium. (iii) Multi-body scatterings may also increase the effective thermal annihilation cross section [64][65][66][67][68]. (iv) A non-trivial thermal evolution in the early universe may depopulate an initially over abundant DM relic density [69].

Non-thermal production
The DM abundance may also be produced non-thermally. If DM is only very weakly coupled to the SM thermal bath and it has not been produced during reheating, DM may be slowly produced via the freeze-in mechanism [39,40]. Ref. [30] discussed the phenomenology of the freeze-in mechanism in the scotogenic model [9] with fermionic DM, one of the examples where DM-nucleus scattering occurs at one-loop level.

Conclusions
Direct detection of DM may not have been observed yet because it is absent at tree level, occurring only at the loop level. In this work we have studied the case of a fermionic singlet DM ψ, which is a simple scenario where DD is naturally induced at one-loop order. The type of scenario considered appears in supersymmetric extensions where the neutralino is pure bino [15] (notice that in this case its mass is typically very heavy, larger than 2 TeV), and also in connection to neutrino masses, in particular in the seesaw model [13] and in some radiative neutrino mass models [8,9,12,80]. We have considered a simplified scenario with a dark sector made of a vector-like (or a SM) fermion and a (complex) scalar. We presented general analytical expressions for the different contributions as well as current limits on the dark sector parameters. We have outlined the possible UV completions of the corresponding penguin diagrams, also those involving SM fields, and we summarize the different possibilities in the following: (i) If the fermion is a SM lepton and thus leptophilic, the DM interactions are generically flavored [5] and there is an interesting phenomenology. There may be new contributions to the anomalous magnetic moment, but the limit is very weak. If there are couplings to at least two different flavors, there are strong limits from LFV, especially for couplings to both electrons and muons. In this case the limits from LFV processes such as µ → eγ and µ → 3e are much stronger than DD. In the absence of one of these couplings DD limits are stronger above a certain DM mass given by the experimental energy thresholds of the DD experiments. In some cases the same particles entering in the DD loop may naturally violate lepton number (specially if the DM couples to the left-handed lepton doublets) and give rise to radiative neutrino mass models such as the scotogenic model with Majorana DM [9] or the generalized scotogenic model with Dirac DM [80].
(ii) If the dark fermion is a right-handed neutrino, it may be a Majorana fermion and an active Majorana neutrino mass term is generated via the seesaw mechanism [81]. As the particles in the loop are neutral, DD is generated via Z and Higgs penguin diagrams [13], which are very suppressed. Although the DM may be assigned lepton flavor and lepton number, there are no strong limits from LFV or lepton number violation beyond those already present in seesaw scenarios. This scenario is normally referred to as the neutrino-portal to DM [13,41,42].
(iii) If the scalar is the SM Higgs, there is mixing between the DM and the neutral component of the fermion in the loop, which generates tree-level contributions mediated by the Z boson and the Higgs. The Z-mediated tree-level DD is expected to dominate with respect to the dipole moment contributions arising at loop level. In fact, elastic Z-mediated contributions are already ruled-out by DD experiments.
While the correct relic abundance is easily achieved in models with DM couplings to SM leptons (or not too heavy right-handed neutrinos), it requires further model-building in the case of DM couplings to vector-like fermions. We have also found that the invisible loop-induced Z and Higgs boson decays may sometimes impose restrictions in the case of light DM.
In this work we studied the prototypical case of fermion singlet DM with the simplest dark sector, where the loop suppression still allows reasonably large DM interactions. Hopefully a positive DD signal in the next years will serve as a motivation and guidance to continue exploring the WIMP DD theory space and its interplay with other beyond the Standard Model probes.

A Larger dark matter groups
In the main part of the text we restricted ourselves to a global U(1) symmetry for a Dirac DM and to a discrete Z 2 symmetry for a Majorana DM. Our results can be easily generalized if the DM forms a larger non-trivial representation of the dark symmetry group and there are multiple degenerate components of the DM multiplet. As the dark symmetry commutes with the SM gauge group it simply leads to an overall factor of to the Wilson coefficients of a DM particle-nucleus scattering, ψ α N → ψ α N , where the Clebsch-Gordan coefficients C γ βα are defined such that the scalar and the two fermions are invariant under the dark sector symmetry: Thus for a general DM candidate with N components ψ α the DD cross section is obtained by summing over the final states and averaging over the initial state and thus Note that a larger dark sector symmetry may lead to multiple DM candidates, which requires to go beyond the discussed scenario, see for instance Ref. [83].
B Direct detection differential cross section and event rate The differential cross section for fermionic DM may be written in terms of NR operators at the nucleon level [52] dσ with the nucleus mass m A and spin J A . The coefficients R X are given in terms of the NR Wilson coefficients c 0,1 i = (c p i ± c n i )/2 and W X denote the nuclear response functions. The explicit forms of R X and W X are given in Ref. [50]. For | q| → 0, the long wavelength limit, W M (0) ∝ A 2 counts the number of nucleons in the nucleus, W Σ and W Σ measure the nucleon spin content of the nucleus, W ∆ measures the nucleon angular momentum and W ∆Σ the interference.
In the literature it is also common to show the differential cross section as the sum of different dipole and charge contributions. Neglecting the Z contributions to SD interactions, which are suppressed with respect to the long-range interactions, and taking d ψ = 0, the differential cross section can be written as [84]: where µ ψA = m ψ m A /(m ψ +m A ) is the DM-nucleus reduced mass and A eff encodes the DM-nucleus couplings (see e.g. Ref. [10]): The first line in Eq. (51) corresponds to the dipole-charge (D-C), the second line to the dipoledipole (D-D) and the third line to the charge-charge (C-C) interaction. F SI (E R ) and F SD (E R ) are the nuclear form factors. c N SI with N = n, p are the relativistic Wilson coefficients at the nucleon level for the operators The vector operator O N,V SI is induced by both interactions with a photon and a Z boson. Once the differential cross section is computed via Eq. (50), the differential event rate per unit detector mass (for a detector with just one type of nucleus A) is given by: where ρ ψ is the local WIMP density, f det ( v) is the WIMP velocity distribution in the detector rest frame and v min is the minimum WIMP velocity required to produce a recoil with energy The velocity distribution in the detector rest frame is related to the velocity distribution in the galaxy frame f gal ( v, t) by a simple Galilean transformation, is the velocity of the Earth in the galactic frame. In our analysis we use LikeDM and refer to [62,63] for the technical details of the different detectors and astrophysical assumptions.

C Expressions for Z and Higgs boson decays into dark matter
The relevant interactions of the DM ψ with the Higgs and the Z boson can be parameterized as 17 and where p µ 2 is the 4-momentum of the outgoing DM ψ. We define x h ≡ m ψ /m h and x Z ≡ m ψ /m Z . The partial Higgs decay width into the DM ψ is non-zero for m ψ < m h /2 and reads: Similarly the partial width of the Z is given by for m ψ < m Z /2. S is the symmetry factor, equal to 1/2 for identical final states (Majorana DM), and equal to 1 for Dirac DM. The coefficients relevant for the decays of the Higgs boson to Dirac DM can be expressed in terms of the Passarino-Veltman functions 17 In the case of radiative neutrino mass models such as the scotogenic model [9] and its variants [38,80], there are extra (lepton number conserving) invisible Higgs boson decays into neutrinos at one loop, which are not suppressed by phase space and could therefore be larger than those into DM.
The mass insertions, m ψ and/or m F are needed in order to flip chirality. We do not report the expressions for the decays of the Z boson, as they are very long and not illustrative. For Majorana DM c V = d V = d A = 0 and the remaining non-zero Wilson coefficients are a factor of two larger,

D Lepton flavor violation and anomalous dipole moments
If the DM couples to SM leptons there may be LFV processes and anomalous electric and magnetic dipole moments. We provide the relevant expressions for DM coupling to either the left-handed SM doublets or the right-handed SM singlets. The results are identical for Dirac or Majorana DM.

D.1 Left-handed lepton doublet
The relevant interaction term for LFV processes is with the charged scalars: The most general amplitude for the electromagnetic charged lepton flavor transition α (p) → β (k) γ * (q) can then be parameterized as [85] A where e > 0 is the proton electric charge, p (k) is the momentum of the initial (final) charged lepton α ( β ), and q = p − k is the momentum of the photon. As is well known, the charged lepton radiative decays are mediated by the electromagnetic dipole transitions in Eq. (63) and the corresponding branching ratio (Br) for α → β γ is given by where with For trilepton decays we consider only the contributions from the photon penguin and from boxtype diagrams, as the Z penguin is suppressed by charged lepton masses. Box diagrams may be the dominant contribution in absence of the contributions from photon and Z penguins. The amplitude from the box diagrams is given by A BOX = e 2 B u(k 1 ) γ α P L u(p) u(k 3 ) γ α P L v(k 2 ) .
The contribution from box diagrams B for − α → − and for − α → − γ − γ + β it is given by with All the external momenta and masses have been neglected. Of course for α → β β β both Eq. (73) and Eq. (74) agree with γ = β.
For µ − e conversion in nuclei we only consider coherent scattering via photon contributions, but include both short-and long-range contributions [86]: 18 L int = − e 2 m µ A L 2 e σ µν P L µ F µν + m µ A R 2 e σ µν P R µ F µν + h.c.
− q=u,d,s g γ LV (q) e γ α P L µ + g RV (q) e γ α P R µ qγ α q + h.c. .  Table 7: The overlap integrals in the units of m 5/2 µ and the total capture rates for different nuclei [86]. The total capture rates are taken from Tab. 8 in [86]. The overlap integrals of 197 79 Au as well as 27 13 Al are taken from Tab. 2 and 48 22 Ti are taken from Tab. 4 of Ref. [86].
The µ − e conversion rate is where the effective vector couplingsg (p,n) L/RV for the proton and the neutron arẽ with g γ LV (q) = e 2 Q q A L 1 .
The coefficients A L, R 1 are given in Eq. (71), and Q q is the quark electric charge of the quark q in units of e > 0. The numerical values of the overlap integrals D and V (p,n) and the total capture rate for each nucleus are reported in Tab. 7 for three different nuclei. As we only consider the photon contribution and thus only couplings to the electric charge of the quarks, there is no effective coupling to neutrons.
Even if lepton flavor is conserved there are processes that can bound the DM interactions with the leptons. Electric dipole moments for the leptons occur in these simplified models only at the two-loop level. However leptonic magnetic dipole moments occur at one-loop order via photon penguin diagrams, similarly to µ → eγ transitions. They receive two independent contributions from the charged scalars running in the loop, which are given by [87]: A R 2 is the diagonal part (α = β ≡ ) of the coefficient given in Eq. (65) and the loop function is defined in Eq. (66). Our expression agrees with Ref. [4]. In the case of the muon magnetic dipole moment, the discrepancy with the SM has the opposite sign and hence the model cannot explain it. However this can be used to (very weakly) bound the model. Electron and tau AMMs do not lead to any relevant constraints.

D.2 Right-handed charged lepton
The relevant interaction term for LFV processes is with the charged scalars: All the expressions are the same as for the left-handed lepton doublets after substituting the right-handed superscript by the left-handed one, i.e., A R 1, 2 ↔ A L 1, 2 , g γ LV ↔ g γ RV and the Yukawa couplings y 1 ↔ y 2 .

E Computation of the relic abundance
In this appendix we review the computation of the relic abundance, see for instance Refs. [20,88]. We use the instantaneous freeze-out approximation which is sufficient for our purposes. The final DM abundance is determined by with λ = [xs σv /H] x=1 and n = 0 (1) for s-wave (p-wave) DM annihilation. The entropy density is denoted by s, with today's value given in terms of the CMB temperature T γ,0 = 2.73 K as s 0 = 2π 2 /45 (43/11) T 3 γ,0 , where we have used N eff = 3. Equating the interaction rate Γ ann for the process ψψ ↔ α¯ β with the Hubble rate, H(T f ) = Γ ann (T f ) we obtain a condition for the freeze-out temperature where m P is the Planck mass, g * is the number of relativistic degrees of freedom at freeze-out (g * = 106.75 in the SM), and g ψ is the DM number of degrees of freedom, which is equal to 2 (4) for Majorana (Dirac) DM. The annihilation cross section may implicitly depend on the freeze-out temperature, and it is useful to factorize out this dependence. Then Eq. (83) can be written in terms of λ as Solving for λ in Eq. (84), plugging it in Eq. (82), and imposing that relic abundance matches the observed value Ω ψ h 2 = 0.12 [89], one can numerically obtain the value of x f . We get values of 23 x f 30 for 10 GeV m ψ 10 4 GeV, which turn out to be identical for Dirac and for Majorana DM. We also note that x f increases roughly logarithmically with the DM mass m ψ .
For a given (m ψ , x f ) pair Eq. (84) allows one to compute the annihilation cross section averaged over velocity, σv , which depends exponentially on x f . For the range of DM masses given above the dependence on the DM mass is very mild. We obtain that the required thermally averaged annihilation cross sections to reproduce the observed DM abundance are in the range 1.8 10 26 σv | D (cm 3 s −1 ) 2.4 and 4 10 24 σv | M (cm 3 s −1 ) 9 for Dirac and Majorana DM, respectively.

F Matching onto non-relativistic operators
We use DirectDM [48] which follows the normalization of the NR operators in Ref. [50] to match our one-loop calculation of DM scattering off quarks onto the NR operators using 3 flavor QCD without running, i.e. the matching occurs at µ = 2 GeV. This is justified as the relevant relativistic operators are renormalization group invariant under one-loop QCD corrections. There are no additional significant contributions, because the particles in the loop are color singlets. There can and neutrons c n 1 = 0.072c d SS − 0.12563c g + 0.0826c s SS + 0.03c u SS (99) c n 4 = −7.176c d AA + 0.248c s AA + 3.008c u