Dark Matter EFT, the Third -- Neutrino WIMPs

Sterile neutrinos coupling only to third-generation fermions are attractive dark matter candidates. In an EFT approach we demonstrate that they can generate the observed relic density without violating constraints from direct and indirect detection. Such an EFT represents a gauged third-generation $B-L$ model, as well as third-generation scalar or vector leptoquarks. LHC measurements naturally distinguish the different underlying models.


Introduction
What particles form dark matter (DM), how do they couple to the Standard Model, and how did they get produced? These questions require a fundamental interpretation framework for a stable new particle from a new physics sector. After defining its hypothetical quantum numbers we can add the dark matter particle to the Lagrangian of the Standard Model (SM) and obtain a perfectly consistent theory. If this Lagrangian is perturbatively renormalizable, the couplings to the Standard Model are given by the quantum numbers of the dark matter agent. If the Lagrangian includes higher-dimensional operators, the interaction through a decoupled new mediator defines a dark matter effective theory. The only constraint we have on these interactions is the measured relic density, combined with an assumed production mechanism in the thermal history of the Universe.
Of many candidates and scenarios, weakly interacting massive particles (WIMPs) [1][2][3][4][5] with their defining thermal freeze-out production still stand out [6,7]. The original papers all use sterile neutrinos as actual WIMPs, and we will return here to this notion. Direct detection experiments are slowly covering the relevant WIMP parameter space [8,9]. To obtain the observed relic density from freeze-out production, the dark matter annihilation rate to SM-particles has to be large. Hence, SM-mediators like the Z-boson or even the Higgs boson are ruled out over almost the entire parameter space. New mediators can avoid these constraints, and if the annihilation process is not enhanced by an s-channel funnel it can be described by a dark matter EFT (DMEFT) [10,11]. The effective Lagrangian can be informally written as L = C/Λ 2 f 2 DM 2 , where f is a SM fermion and DM the DM particle. We can then translate the observed relic density Ω DM h 2 ≈ 0.12 into a typical dark matter annihilation cross section [12] σ ann v rel ! ≈ 1 600 TeV 2 ∼ The relation between the dark matter mass and the new physics scale becomes Λ 2 C m DM = 6.9 TeV This indicates that large mass ratios R = Λ/m DM require light new physics to produce the correct relic densities, while large Wilson coefficients C can alleviate this pressure. Too small values of R will cast doubt on the EFT assumption, but we can assume a typical range R = Λ m DM = 3 ... 10 and C 1.
This picture of the DMEFT is consistent and describes not only the relic density, but also direct and indirect detection.
Matching to the full theory we trade the new physics scale Λ and the Wilson coefficient C for mediator mass and coupling, Λ ≈ m med and C 2 ≈ g 4 . This matching allows us to include collider measurements in a global analysis, even if these measurements require an on-shell production of the mediator. At the same time, matching to a full theory brings in constraints from the construction or the consistency of the model. For instance one needs to address flavor observables or assure perturbatively small values of g whenever the model should be valid to high scales. If there exists no UV-complete theory completing a consistent EFT approach, such an EFT is utterly pointless. From this perspective, a fully consistent EFT approach is less trivial than one may think [13].
In this paper, we build on the recent appearance of 4-fermion effective interactions in the neutrino sector [14][15][16][17][18], beyond oscillation analyses [19][20][21][22] and consider a DMEFT approach to a sterile neutrino. After the original DMEFT papers [10,11] and a largely unsuccessful attempt to turn DMEFT into a global fitting framework [13,23], this third attempt tries to accommodate these earlier complications. We let it couple only to the third fermion generation, which avoids many of the arguments that make consistent DMEFTs difficult, and thus allows for a consistent EFT approach for a variety of complete models. First, considering the relevant higher-dimensional operators, we demonstrate that the relic density via freeze-out can be generated, while at the same time direct and indirect detection limits can be obeyed without loosing the validity of the EFT approach. Next, we investigate possible UV-completions of this framework: • the anomaly-free gauged Abelian (B − L) 3 symmetry. The dark matter sterile neutrino couples to a Z which in turn couples to third generation fermions; • a scalar third generation leptoquark, which couples to the dark matter sterile neutrino and right-handed top quarks; • a vector third generation leptoquark, which couples to the dark matter sterile neutrino and right-handed bottom quarks.
These models face different constraints in particular from LHC searches, as well as direct and indirect detection.
The paper is built as follows: In Section 2 we discuss the EFT of fermion singlets, which includes a variety of higher-dimensional operators. Focusing on couplings to third generation particles, we identify the operators relevant for DM freeze-out, as well as direct and indirect detection. We shown that our EFT scenarios allow for the successful DM generation in agreement with all constraints. In Section 3 we discuss some issues with consistent EFT scenarios. We then match our EFTs with explicit models, specifically the gauged Abelian (B − L) 3 as well as scalar and vector leptoquarks. The models are confronted with additional constraints from colliders and flavor, before we conclude in Section 4.

νDMEFT framework
For our νDMEFT we assume a sterile neutrino N of mass m N 10 GeV as the DM agent. * Its mixing with active neutrinos must be very small or forbidden by symmetry, for instance if N is the only particle odd under a Z 2 symmetry. This forbids in particular a Yukawa coupling with the SM Higgs, which would make it unstable [24,25]. The standard seesaw mechanism as one possible origin of Majorana neutrinos N implies that other sterile neutrinos would have to be heavier, so they decouple and decay before N freezes out. For instance, two heavy sterile Majorana neutrinos are sufficient for the seesaw mechanism to generate light neutrino masses (one of which would be massless), and leave a super GeVscale sterile neutrino as a DM candidate. However, this realization is only one possibility, there are countless ways to generate active neutrino masses. We will also consider the Dirac option for N in the course of the paper, and will comment on the difference to the Majorana case when appropriate.

Operator basis
We describe a sterile Majorana neutrino as a 4-component spinor where N R and N c R are 4-component chirality eigenstates. Assuming only SM gauge symmetries, the renormalizable operators are given by the SM-Lagrangian plus a sterile neutrino contribution, We extend this Lagrangian by D5 and D6 terms, At dimension five the SMEFT [26] only features the Weinberg operator, while the neutrino extension allows for two additional terms [27], Here B µν is the hypercharge field strength. For only one Majorana neutrino, the magneticdipole-like O  generates a mass correction for the sterile neutrino, but also Higgs (portal) couplings of the vertex forms N N h and N N hh. Whenever a pair of sterile neutrinos may annihilate into a SM fermion pair, this operator is also produced at one-loop level. This is the case for all models studied in this work. Therefore, we will occasionally return only to O (5) N H out of the D5 operators in the following.
An operator basis at dimension six, including sterile Majorana neutrinos, is discussed in Refs. [28] and [29], see also Ref. [21]. We reproduce the Warsaw basis in Table 1 for the 4-fermion operators and Table 2 for mixed fermion-boson operators. In the upper sections of the tables we list the SMEFT operators, while in the lower sections we add the operators involving sterile neutrinos (νSMEFT). We leave out operators not connected to fermions, as well as lepton or baryon number violating operators unless originating from N , and the operator O N 4 = N c R N R N c R N R which violates lepton number by 4 units. We note that the first sterile-neutrino-gluon operator arises at dimension seven, where one or both gluon field strengths can be exchanged by their dual.
In the Dirac case, all operators violating lepton number are absent, and in Eq.(5) we need to replace the Majorana kinetic and mass terms by their Dirac fermion counterparts. The operator O as will be explained with the direct detection constraints. Apart from that, we consider the same D6 operators coupling only to the right-handed component of the sterile neutrino, although in principle interactions also with the left-handed components could be present in this case. Other differences arise when annihilation cross sections are important.

Constraints
We need to gather the existing constraints for our νDMEFT, assuming that N is the only field odd under a Z 2 symmetry. In this section we will use the pure EFT approach and  derive limits on the Wilson coefficients and cut-off scale when N couples to either τ R , t R , or b R . For each of these three cases we consider the operator O N τ (t,b) ≡ O 33 N e(u,d) with Wilson coefficient C N τ (t,b) . Relevant constraints come from the relic abundance, as well as direct and indirect detection. The results are shown in Figure 1 for the Majorana and Dirac cases.
Relic abundance For a valid DM candidate we need to ensure that Ω N h 2 ≤ Ω DM h 2 = 0.12. The key ingredient is the annihilation process where the relevant particles in the final state are defined by the mediator or the corresponding effective operator. In all our cases the annihilation goes into fermion-antifermion pairs. The EFT framework with Λ m N naturally matches the non-relativistic freezeout scenario, provided the 2 → 2 annihilation rate is large enough to predict the observed relic density. We check the correct relic density using micrOMEGAs [30][31][32][33]. The black lines in Figure 1 correspond to Ω N h 2 = 0.12, while the gray area overshoots the observed DM density. Note the different position of the DM-line for Dirac and Majorana fermions, reflecting the different annihilation cross sections. When kinematically allowed, we find the cross sections for the operators O N f with a right-handed singlet f = τ, t, b to read for Majorana DM and up to first order in average squared relative DM velocity v 2 rel by direct computation. For Dirac DM the annihilation rate is Here we include the color factor N c = 1 for f = τ and N c = 3 for f = t, b. In the Dirac case, we drop terms of order For the lepton case with operator O N we set N c = 1 and replace the masses m b with m τ and m t with m ντ = 0.
Direct detection The most promising channel to directly detect multi-GeV dark matter such as the traditional WIMP is elastic scattering off nuclei, searched for at experiments such as XENON1T [34]. To contribute to this scattering, the new physics must couple DM to light quarks or gluons at nuclear energy scales. This may happen in various ways at tree level or loop level. Mapping these interactions on non-relativistic nuclear scattering theory allows one to compare the predicted scattering cross section with experimental bounds, see e.g. Ref. [35] for a description of a mapping from the UV theory to the nucleon-level theory.
For our νDMEFT the scattering off light quarks is described by four-fermion operators O N f , just as in any WIMP dark matter scenario with an underlying Z 2 symmetry [36]. For this to be detectable, the operators O N q , O N u , or O N d with quark flavor indices 1 or 2 must be present. However, even if their Wilson coefficients vanish at the weak scale, there will be a non-vanishing coupling at the nuclear scale induced by RG-running. Technically, this means that one needs to map the Wilson coefficients to the appropriate EFT of the SM extended by a right-handed neutrino below the weak scale. This mapping has been discussed in Ref. [21] for neutrino interactions and recently been given completely [37]. We use runDM [35] to check whether the running-induced N -nucleus coupling is or is not within the reach of current experiments. We find that in the Majorana case, when the N -bilinear reduces to an axial current, N γ µ γ 5 N , only very weakly constrained operators of non-relativistic scattering theory are generated. We typically find couplings of order 10 −4 to 10 −9 for operators 5, 9, and 10 defined in Ref. [38], well below current XENON limits [39]. In the Dirac case, however, the N -bilinear involves a vector coupling 1/2N γ µ N . In this case, the more stringently constrained O 1 in non-relativistic scattering theory is generated, where in the scenarios of τ -and b-coupling the proton coupling dominates and in the scenario of t-coupling the neutron coupling dominates. To apply bounds on this operator from XENON100, we use the Mathematica package provided in Ref. [38] to calculate the event rates in the detector as a function of the non-relativistic coefficient. This coefficient can be mapped to the new-physics scale Λ with runDM. We then use these rates along with the detector response functions provided by XENON [39] to calculate the expected event numbers. Since no excess over background was observed, we use this to obtain 90% CL bounds via a simple chi-square function. According to Ref. [39], the detector response functions we used should provide conservative but not-too-far-off limits. We have confirmed this by reproducing the limit for O 1 in Ref. [39] which in our case is about an order of magnitude more conservative. Only in the t case, this procedure yields a lower bound on the mass of Dirac DM when combined with the relic density constraint, as can be seen from Figure 1, namely m N ≥ 1549 GeV.
Two more operators may contribute to a light-quark coupling, namely O HN and O (5) N H . The former will not be encountered in this paper. The latter introduces a Higgs portal coupling N N h, as mentioned after Eq. (7). This Higgs portal coupling is severely constrained by direct detection. To quantify this, we calculate the spin-independent elastic N -nucleus scattering cross section in the zero-momentum-transfer limit [40] with the reduced nucleon mass µ x = m N m x /(m N + m x ), and for x = p, n, where we take for f x T q and f x T G the same values as applied in Ref. [41]. We compare the prediction for the Higgs-mediated direct detection cross section with the latest XENON1T results [34] to derive 90% CL lower limits on Λ/C (5) N H . In the sensitive range of DM masses (6-1000 GeV), we find Λ/C (5) N H > 6 · 10 4 GeV to be consistent with all masses and Λ/C (5) N H < 2 · 10 3 GeV to be excluded for all masses. In Figure 2 we show the mass-dependent limit on Λ/C (5) N H derived under the assumption that the entire cross section is determined by O N H is so strongly constrained, that it is often justified to neglect the annihilation channel N N → h when considering a freeze-out through four-fermion interactions. Strictly speaking, the Higgs channel may still play a role if the DM mass is close to the resonance m N ∼ m h /2. We do not consider this possibility further in this work. On the other hand, this strong constraint allows to test the models presented in Section 3 in direct detection even when the operator is only generated at loop level. In the EFT scenarios considered in Figure 1, we assume Λ/C (5) N H > 6 · 10 4 GeV such that there is no restriction in the parameter space from direct detection. In the Dirac case, we included a factor of two in Eq. (9), to keep the direct detection cross section equal to the Majorana case for the same coefficient C N H introduced in Eq. (7). The area is excluded at 90% CL.
In addition to quarks, DM may interact with gluons in the nucleus. This interaction is described by operators like O (7) gN of at least dimension seven. We ignore it for the D6-Lagrangian considered here. If DM couples to quarks, this operator is necessarily generated via loops. It can be more relevant than the lower-dimensional four-fermion operators, for instance for the Higgs portal, where the light quark Yukawas are tiny and the top loop does not decouple. For other mediators coupling only to third-generation fermions, the power suppression from the loop will weaken the direct detection limits. As an example, for the model described in Section 3.2 a Higgs portal and a gluon coupling, both generated at one loop, are the relevant contributors to direct detection.

Indirect detection
The νDMEFT operators relevant for indirect detection are typically the same as the operators leading to the observed relic density. A significant mismatch between these two processes occurs only when the annihilation process involves a narrow s-channel resonance, which cannot happen in the EFT description.
We consider the search for gamma rays from DM annihilation in dwarf spheroidal galaxies (dSphs) by Fermi-LAT [42]. The expected gamma ray flux in the energy range between E min and E max , and allowing for more than one annihilation process (j), is given by where σv rel j denote the velocity-averaged annihilation cross sections and dN γ,j /dE γ the corresponding photon spectra. In addition, this formula includes a J-factor for each considered object. The Fermi-LAT and DES collaborations published bin-by-bin likelihoods for the gamma ray spectra of a number of dSphs [42]. We compare the predictions from Eq.(16) with these bin-by-bin likelihoods for each dSph to derive limits on the annihilation rates. The inputs include a set of dSphs, their J-factors with corresponding uncertainty, and dN γ,j /dE γ . The calculation of the photon spectra requires particle showers and we employ the spectra provided in Ref. [43].
We include all 19 dSphs with kinematically determined J-factors in Table 1 of Ref. [42]. For the J-factors we employ the more recent results in Ref. [44] and estimate the uncertainties by the 68% confidence interval given in that reference. For a given m N we then use the best of the 19 limits on σv rel j using the upper 68% CL bound on the J-factors as the best-case limit, and the best out of the 19 limits using the lower 68% CL bound as the worst-case limit. This defines an exclusion band between those two values. We also note that the J-factors from Ref. [44] are in most cases slightly lower than previous estimates, so our limits are slightly weakened. The limits obtained by us with this procedure are shown in Figure 1. The light (dark) blue areas correspond to the best-case (worst-case) scenario of the J-factors being at the upper (lower) boundary of their 68% CL region. Note the different position of the indirect detection areas for Dirac and Majorana fermions, which is caused by the different form of the annihilation cross sections, see above. From the intersection of the black line with the blue areas, we can derive the following lower limits on the DM mass in the Majorana case m N ≥ 6 -14 GeV for bb , 6 -10 GeV for τ τ , and in the Dirac case In the tt case, there is no intersection and the whole parameter range is consistent with indirect detection.
Even though we could strictly truncate the EFT expansion at dimension six, we add another indirect detection channel for the Majorana case, namely the aforementioned D7 operator O (7) gN . This is motivated by existing literature showing that a simple model generating the tt and bb scenarios features a significant contribution to the gamma ray spectrum from loop-induced N N → gg annihilation [45]. It is significant only because the Majorana annihilation cross section into fermion pairs is suppressed at present times, since the second term in Eq.(11) involves a factor v 2 rel ∼ 10 −6 . We can apply the results from Refs. [46,47], which give expressions for the loop-induced annihilation rates of neutralinos into two gluons, mediated by tops and stops or bottoms and sbottoms, respectively. The stop corresponds to the scalar leptoquark of our UV-complete theory in Section 3.2, which generates the EFTtt case discussed here. Analogous results for vector leptoquarks, corresponding to thebb case as outlined in Section 3.3, are not available in the literature. Following the same procedure as above and setting v rel = 10 −3 , we include this channel in the derivation of indirect detection constraints. The constraints for a combination of both fermionic and gluonic annihilation channels are seen as the orange continuations in the bb and tt plots of Figure 1. The effect is less pronounced in the tt case, since the first term in in Eq. (11) is not as suppressed for top quarks with their large mass. In both cases the effect is not strong enough to affect the lower limit on m N from indirect searches discussed above. In the τ τ case, this channel is expected to be irrelevant due to τ leptons being color singlets.

νDMEFT representing models
While in some instances effective theories can safely be viewed as stand-alone theories, we know that for a DMEFT it is crucial that the effective theory can also be justified as a low-energy approximation to classes of UV-complete models [13,48]. Otherwise it would be an utterly pointless analysis. This means we need to compare relevant theory predictions between the DMEFT and corresponding models with propagating mediators. The observables we need to consider are DM annihilation predicting the freeze-out relic density, direct detection, indirect detection, flavor constraints, and collider searches. Some scenarios where the LHC has obvious potential do not provide enthusiastic support for global DMEFT analyses [13,23]: • tree-level s-channel vector mediator coupling to first-generation quarks: the mediator is strongly constrained by LHC resonance searches. In the allowed parameter range the observed relic density requires an s-channel funnel and invalidates the EFT approach; • loop-level s-channel scalar mediator coupling to third-generation quarks at tree level and to gluons at loop level. DM can annihilate into heavy quarks, gluons, and mediator pairs, challenging the EFT picture. Collider searches for resonant mediator production lack sensitivity; • tree-level t-channel scalar mediator coupling to first-generation quarks: the mediator is constrained by LHC pair production. In the remaining parameter range the annihilation rate is too small to explain the observed relic density; • loop-level t-channel mediator coupling to third-generation quarks at tree level and mediating a DM coupling to gluons at loop level. DM can annihilate into heavy quarks and into gluons, challenging a fixed-dimension EFT approach. Nonetheless, we revisit this class and identify suitable regions of the parameter space where the D6-EFT is consistent except for the LHC signatures.
For these classes of mediators the main observation is that mostly the LHC constraints on the mediators make it hard to predict the observed relic density with a DMEFT.
The νDMEFT adds a new set of UV-models with the promising feature that they include tree-level mediators which only couple to third-generation fermions. For an EFT approach to make sense, it needs to represent more than one specific UV-model. Examples for such models are a gauge extension only coupling to third-generation fermions in Section 3.1, a scalar third-generation leptoquark coupling to t-quarks in Section 3.2 (which represents the loop-level t-channel mediator class), and a third-generation vector leptoquark coupling to b-quarks in Section 3.3. For these three UV-completions we apply the EFT framework to the relic density, direct and indirect detection. For the LHC we accept that on-shell mediator production cannot be analyzed in the EFT framework [49], so we translate the EFT constraints from cosmological constraints into full-model parameter regions which should be tested at the LHC.

Gauging third-generation (B − L)
For gauge extensions it is attractive to start with global symmetry groups which do not develop anomalies when gauged [50,51]. Because anomalies do not require cancellations between generations we gauge B − L for the third generation only and avoid most constraints while concentrating on flavor physics [52] or dark matter features [53]. If we neglect generation mixing [53] we need one additional SM-singlet scalar Φ, which carries (B − L) 3 charge +2 and breaks the additional U (1) symmetry. Models based on B − L are tailor-made for Majorana fermions, so we do not consider the Dirac option here. With a discrete Z 2 symmetry under which only N R is odd we find the Lagrangian where µ 2 H and µ 2 Φ are chosen positive, f = l τ , τ R , q 3 , t R , b R , and The last term in the first line of Eq. (19) is the only Yukawa coupling involving Φ. The Z 2 and (B − L) 3 symmetries forbid many terms discussed in Section 2.1. We assume for the kinetic mixing parameter that 1 in order to satisfy stringent experimental constraints, as detailed later in this section. The field Φ is a SM singlet, hence only kinetic mixing is present.
Masses and mixing Below the weak scale and in unitary gauge, the Higgs field and the additional U (1)-breaking scalar Φ can be expressed as where v and w denote vacuum expectation values of H and Φ respectively. We omit the procedure to obtain the physical masses and currents here and refer to the general case discussed e.g. in Ref. [51]. If w v we can expand the VEVs in λ HΦ 1, We can describe the small mixing of h and φ with an angle θ ∈ [−π/4, π/4], given by The new physics mass spectrum is then wherem X is the mass of the non-canonically normalized fieldX below the U (1)-breaking scale. In the limitm X v the masses of the physical Z and Z fields in the broken electroweak phase are where m Z,0 = gv/2c W and we drop terms of order 4 and m 4 Z,0 /m 4 X . The limit m N m Z then corresponds to y g X 1, if the U (1) (B−L) 3 is to remain perturbative.
The additional Z-Z mixing induced at the one-loop level is calculable, but in our case it requires a renormalization of . In practice, the renormalized then becomes a free parameter which must be fixed by a measurement and can be constrained from electroweak precision data, depending on m Z [54].
Operator matching at the weak scale We assume that the symmetry breaking scale of the (B − L) 3 is above the weak scale, so in the matching there is no scalar mixing. In Eq.(21) the heavy VEV is now just w = µ Φ / √ λ Φ . The scalar and vector masses are The interactions of φ, excluding self-interactions, then are Any lepton-number violating process, in particular the operators in Eq.(7) must involve an insertion of w.
We start with the D5 operators in Eq.(7). Since ν and N do not mix, the lepton number violation in N is not carried over to the active neutrinos and O νν does not appear. Next, O N B is forbidden for a single Majorana neutrino. However, there exists a tree-level contribution to O (5) N H proportional to yw ∼ m N , This operator has been discussed in Section 2.2 to be potentially relevant for direct detection. However, when the mass of φ is moderately large, this operator is too suppressed to leave a detectable signal, as will be quantified below.
The 4-fermion operators O f f with f, f = l, e, N, q, u, d are generated at tree level by integrating out the heavy vector X. Depending on the generation indices this requires an additional 0, 1, or 2 powers of . The tree-level contribution to the corresponding Wilson coefficient reads where q f X,Y are the charges of fermion f under the U (1) X and hypercharge symmetries. If neither f α nor f β carries a (B − L) 3 charge, i.e. α, β = 1, 2, the contribution is suppressed by a factor of 2 . This leads to a natural hierarchy in contributions to the operators. The operators of type ψ 2 H 3 in Table 2 are generated via the quartic SM-Higgs coupling. An additional contribution is generated at tree level by an internal φ-line, The operators of type ψ 2 H 2 are generated by B and X exchange for singlet bilinears and by W exchange for triplet bilinears like O Hl . Whenever X is the mediator their suppression is at least /m 2 X , since the X-Higgs coupling is generated by kinetic mixing with B. The same is true for O HN . Finally, the L-violating D6 operator O N 4 is also generated by scalar Φ exchange. We will generally assume that the scalar Φ is heavier than N and that the scalar mixing λ HΦ , as well as the Z-Z mixing are small. This leaves us, to 0th order in λ HΦ and , only with D6 four-fermion operators of the third generation. The following Wilson coefficients are non-zero: where LL = ll, N N, ee, N l, le, N e, LQ = lu, ld, N q, N u, N d, qe, eu, ed, and QQ = uu, dd, qu, qd, ud, and flavor indices 3 are understood. The EFT scale is given by Λ = m Z . We will further consider C N H only in the context of direct detection in order to constrain the scalar mixing. For these operators we need to discuss the relevant constraints at EFT level and specific to a given UV-completion.
Relic density (EFT) We implement the model via FeynRules [55] and calculate the DM relic density after freeze-out using micrOMEGAs [30][31][32][33]. In analogy to the single operator scenarios shown in Figure 1, we also calculate the relic density for an EFT defined by the Wilson coefficients in Eq. (30). The result in terms of the mass parameter Λ/ √ c is shown in Figure 3 along with constraints from indirect detection and lepton flavor universality to be discussed below.

Direct detection (EFT)
Looking at direct detection constraints, the D6 operators discussed above are extremely poorly constrained, again, since the running-induced nucleon couplings are small. Direct detection, however, directly probes O N H , which is determined by the combination of heavy scalar mass and scalar mixing, m φ / √ λ HΦ . In Figure 4 we show the XENON1T limits [34]. The advantage of the EFT approach is that we may immediately recast the bounds on the Wilson coefficient shown in Figure 2 onto bounds on the scalar mass and mixing using Eq. (27).

Indirect detection (EFT)
We discussed our method for evaluating indirect detection constraints arising from the four-fermion operators in Section 2.2. The result of this  Figure 1, which provide stronger limits than the τ + τ − case. In the Majorana case, the blue indirect detection limits leads to a lower limit on the DM mass in the range of m N ≥ 7-11 GeV, while in the Dirac case the limit is m N ≥ 11-29 GeV, cf. Eqs. (17) and (18).

Lepton universality (EFT) Since in this EFT new interactions among four SM fermions
arise only for the third generation, flavor observables can be used to constrain the parameter space. In this case we face constraints from lepton flavor universality tests. In particular, the operators O lq , O ld , O qe , O ed coupling only to the third generation will affect b-meson branching ratios. For the Υ(1S) decay into leptonic final states we use the formula [56,57] where α denotes the fine-structure constant, Q b = −1/3 the b-charge, R n (0) the nonrelativistic radial wave function at the origin, and contains the kinematics. The only appearance of the lepton mass is in K . This is why the SM expectation for the ratio R of decay widths is simply Focusing on the third generation, the SM-prediction R τ µ = 0.992 is consistent with the BaBar measurement R τ µ = 1.005 ± 0.013(stat.) ± 0.022(syst.) [57]. Adding statistical and systematic errors quadratically, we use R τ µ = 1.005 ± 0.026 to constrain the νDMEFT.
In this EFT the operators O lq , O ld , O qe , O ed inducing 4-point interactions between two b and two τ all have the same Wilson coefficients and therefore add to a vector-like interaction . Contours of g X producing the correct relic abundance are shown in black, while the dashed lines represent the EFT limit upon identifying Λ = m Z and c = g 2 X with c defined in Eq. (30). For a given value of the purple areas are ruled out at 95% CL. The area with m Z < 3m N is drawn to signify where the EFT language is valid, namely in the upper left corner. Above g √ 4π ≈ 3.5 the theory ceases to be perturbative.
Their contribution to Υ decays adds directly to the photon-mediated contribution and modifies the branching ratio to Since the effective operator is limited to third-generation fermions it predicts We can derive 68% CL limits from the above-mentioned BaBar limit and find Λ/ √ c ≥ 224 GeV since c = g 2 X > 0. We show the limit as the green area in Figure 3. From the intersection of the green boundary with the black relic density line, we deduce a lower limit on the DM mass consistent with LFU constraints in the EFT regime. In the Majorana case, this yields a stronger lower bound on the DM mass than indirect detection, namely m N 17 GeV, however, at lower CL.
LHC and other constraints on the UV-model Turning to the full model, we plot in Figure 5 the m N -m Z plane and fix at each point the gauge coupling g X such that we predict the observed relic density. This way we reduce the parameter space (m N , m Z , g X ) to two dimensions. We compare this in the same plot to the EFT limit by reducing also the EFT parameter space (m N , Λ, c) to the m N -m Z plane by identifying Λ = m Z and fixing c through the relic density. For the validity of the EFT approach we require m Z < 3m N , to ensure that annihilation cannot proceed through an intermediate Breit-Wigner propagator. We note that this definition of the validity of the EFT approach is process-dependent and makes the minimal assumption that a propagating mediator does not contribute for instance to DM annihilation [58].

Using
√ c = g X from Eq.(30), we plot contours of fixed g X as solid lines and contours of fixed √ c of the same values 0.3, 1.0, and 3.5 as dashed lines to illustrate the slight deviations of the EFT limit from the exact model. We find that as long as m Z > 3m N , the freeze-out is nicely described by the effective four-fermion operators. Since for m f m N and v rel ∼ 0.1 the annihilation cross section σv rel M N f discussed in Section 2.2 around Eq.(10) is insensitive to m f , the relative importance of the annihilation channels is simply σv rel τ + τ − ≈ 2 σv rel ντ ντ ≈ 3 σv rel bb ≈ 3 σv rel tt (37) if m N m t , and without the last approximate equality when m b m N m t . Here, the factor of 2 is due to only left-handed tau neutrinos existing, and the factors of 3 are simply a result of the color factors and Wilson coefficients. Around m Z ∼ 2m N , since Z is an s-channel mediator, N N → Z becomes a relevant and immediately resonant annihilation channel during freeze-out. Therefore, we exclude the region m Z < 3m N in our plot, since we investigate only the regions described by the EFT. The reader interested in this region is referred to Ref. [53]. We now explain which generic constraints in the EFT picture discussed above and which model-specific constraints give rise to the excluded regions in Figure 5.
1. Perturbativity: at g X > √ 4π ≈ 3.5, the theory becomes non-perturbative. Hence our simple tree-level matching becomes unreliable. While this is not strictly excluded, the EFT representing fundamental models turns into a collection of operators and our interpretation ceases to be useful in this regime.
2. Collider production: for a significant portion of the parameter space, m Z is small enough to be produced on-shell at the LHC. Therefore, the ATLAS and CMS searches for spin-1 resonances produce relevant constraints. It turns out that the channel Z → τ + τ − is the most sensitive one for models such as (B − L) 3 , where the DM mediator couples predominantly to third generation SM fermions [53,59].
Currently, ATLAS provides the most stringent limits on the pp → Z → τ + τ − cross section [60]. Details on the calculation of this contour are delegated to Appendix A.
3. Direct detection: the discussion is carried over from the EFT case: The runninginduced coupling of DM to light quarks is too small to be relevant, such that the combination of scalar mass and scalar mixing m φ / √ λ HΦ is likely the strongest source of a direct detection signal and can be constrained by XENON1T as shown in Figure 4. Assuming this to be satisfied, there is no exclusion from direct detection shown in Figure 5 6. Z-Z -mixing: as noted earlier, the Z mixing is a free parameter that must be determined by experiment. Any heavy Z that mixes with the Z boson can be  [21]. Our convention is Q = I 3 + Y /2. The cases considered here are S 1 and U 1 , generating only effective interactions of our DM fermion N with right-handed up-or down-quarks.
constrained by its effect on Z observables. Along with Ref. [53], we employ the 95% CL bound on as a function of m Z obtained in Ref. [54]. To illustrate the influence of these constraints, we show in Figure 5 two plots, one with = 0.01 and one with = 0.1. As can be seen from the purple region, for = 0.01 only a narrow band around m Z ∼ m Z is ruled out, leaving open parts of the parameter space below m Z < 200 GeV. The fact that for = 0.1 masses m Z 320 GeV are ruled out completely closes this window for larger mixings. For completeness, we note that when 0.005 all masses m Z are allowed.
Because this will become relevant later, we briefly comment on the hypothetical case where the Z only couples to bottom or top quarks. We still obtain limits, for instance, by taking advantage of ATLAS searches for resonant associated scalar production with a subsequent decay to the same quarks, pp → qqqq where q = b, t [61,62]. Details on this evaluation are delegated to Appendix B. We find a miniscule constraint in the top case, as shown later in Figure 6.

Scalar leptoquarks
For a second UV-completion we turn to leptoquarks [63], when the leptoquark is a scalar or a vector. Leptoquarks have recently re-surfaced to explain flavor anomalies [64,65]. The general Lagrangian including all possible couplings between two fermions (SM and sterile neutrino) and one scalar or vector leptoquark was given in Ref. [21], generalizing the list from Ref. [66]. We can distinguish leptoquarks with fermion number F = 3B + L = 0 and F = 2: Here S 1 and U 1 only couple to the SM if there are sterile neutrinos. We will focus on those two possibilities. We list the quantum numbers of the leptoquarks in Table 3, along with the νSMEFT operators that they induce upon integrating out a virtual leptoquark [21]. We will assume here baryon number violation, which rules out the operators which, together with Eq.(38) may lead to proton decay. Usually, baryon number conservation is assumed to avoid these strong bounds. Finally, we choose to ignore Higgs-portal-like terms of the form XX † HH † , where X is any of the leptoquarks.
In a first step we focus on the scalar leptoquark S ≡ S 1 ∼ (3, 1, −4/3). It is rather secluded from the SM since it has only one interaction term with a SM fermion and sterile neutrinos. It couples only to right-handed up-type quarks, in our case this is the top. Since S then has top-like quantum numbers, it has the features of a top squark in supersymmetry [67]. The DM phenomenology of this framework has been studied in detail in Ref. [68]. It also corresponds to the model considered in Ref. [45], if we fix the involved SM fermions to top quarks.
The Lagrangian relevant for S reads where Besides the coupling to the right-handed up-type quarks, the coupling to gluons is entirely determined by the strong gauge coupling g s [69]. We need to assume that S is odd under the stabilizing Z 2 , because otherwise the Yukawa term in Eq.(41) is not allowed. This type of interaction is possible both for Majorana and Dirac sterile neutrinos and in both cases only the right-handed components will interact.
Operator matching Straightforwardly, for processes at momentum p, for m S p the effective point interaction after integrating out S reads By Fierz, this can be rephrased as the operator O N t ≡ C 33 N u with when we identify Λ = m S . Hence, the effective interaction mimics a vector-mediated tt ↔ N N interaction and we can compare the scenario to a Z coupling only to t R and N R inducing the same operator. At one-loop level, both O N H and O gN are generated. We only state here the formula for the effective Higgs coupling calculated in Ref. [70] for the Majorana case in the terminology of Ref. [68] (the relation reads g hχχ /v = 2C (5) N H /Λ) and in the limit of √ s, m N m S , where y t is the SM top Yukawa coupling, r = m 2 S /m 2 t and the functions F (r), G(r), and H(r) are given in Eq.(10) of Ref. [68]. This way the effective Higgs portal coupling is not an independent parameter, but determined by m N , m t , and m S setting the freezeout conditions. We omit the explicit Wilson coefficient of O (7) gN , since it is only relevant for indirect detection in our EFT region and in that case we may use the explicit model formulae from Refs. [46,47]. Let us discuss the constraints arising in this scenario, and compare to the EFT discussion in Section 2.2.
Relic density We implement the model via FeynRules [55] and calculate the DM relic density after freeze-out using micrOMEGAs [30][31][32][33]. In Figure 6 we plot the m N -m S plane and fix at each point |x t | such that the observed DM relic density is obtained. We compare this to the EFT limit of the model described by the matching Eq.(44) together with Λ = m S . This means that for each contour of constant |x t | we associate a dashed EFT contour with constant |C N t | = |x t |/ √ 2, where the m S coordinates are given by and the fraction of the right-hand side corresponds to the black curves in third row of plots of Figure 1 corresponding to the Majorana and Dirac cases. In this case we find that as long as m S > m N , the freeze-out is nicely described by the effective four-fermion operators. This is again in contrast to the Z , where already close to m Z = 2m N the EFT limit is broken. Since S is a t-channel mediator, the secluded DM annihilation into a pair of mediators is possible outside the EFT regime.
LHC and other constraints on the UV-model We now explain which generic constraints in the EFT picture discussed in Section 2.2 and which model-specific constraints give rise to the excluded regions in Figure 6.
1. Perturbativity: again, we indicate in gray the region when the coupling becomes |x t | > √ 4π ≈ 3.5, where the theory becomes non-perturbative.
2. Collider production: collider signatures of leptoquarks are classified in Ref. [71]. The main production channel of our leptoquark is gg → SS † and does not depend on the coupling x t , but is entirely determined by the strong gauge coupling g s and the leptoquark mass m S . The value of the coupling x t is only relevant in so far as it is assumed that the decay width is sufficiently small to consider the leptoquarks to be produced on-shell. The most stringent bound on m S can be derived from a CMS stop search [72] (see Ref. [73] for a very similar ATLAS search). We use the analysis that assumes the stop decays into a neutral neutralino and a top, which essentially is our signature. In the case of m N 400 GeV this excludes m S 1200 GeV at 95% CL. We use the exact exclusion contours in Figure 6, they are shown in red.
3. Indirect detection: as can be concluded from the third row of Figure 1 (blue regions), current constraints from annihilation in dSphs are consistent with the full DM mass range that we consider, if only annihilation into fermion pairs is considered. From Ref. [45], however, we know that annihilations into gluon pairs are also relevant in the Majorana case for certain parameter configurations. As described in Section 2.2, the N N → gg annihilation channel becomes relevant at the edge of our DM mass region near m N ∼ 10 4 GeV [68]. The additional contribution from gluons to the expected gamma ray signal is, however, not sufficient to probe the interesting region where Ω N h 2 = 0.12, as can be seen from the distance between the orange region and the black line in the bottom left plot of Figure 1.
gN . In the EFT region m S 3m N , the Higgs-coupling dominates [68]. Therefore, we again apply the more recent XENON1T bound on the effective Higgs operator shown in Figure 2 using the expression for the one-loop Higgs coupling Eq. (45). We find that the limits are not strong enough to probe the parameter space with m N ≥ 160 GeV shown in Figure 6. The bound on the WIMP-nucleon cross section would need to be improved by a factor of about 100 to cut into our parameter space, and a factor 10 4 improvement on the WIMP-nucleon cross section would be required to rule out all of the mass range m N < 1000 GeV. This means that we cannot expect XENONnT to produce notable constraints on the EFT region of this model [74].
To briefly summarize, there remains parameter space for a neutral fermion singlet that couples via scalar leptoquarks to t-quarks, consistent with an EFT description. While the leptoquark mass should exceed the TeV-scale, neutrino masses can be as low as m t , in the Majorana case. In the Dirac case direct detection limits exclude DM masses up to 1549 GeV. Compared to a Z which couples only to t R and generates the same four-fermion EFT operator, we observe that the collider constraints are stronger in the leptoquark case.

Vector leptoquarks
As a third UV-completion we consider the vector leptoquark U ≡ U 1 ∼ (3, 1, 2/3). It is also secluded from the SM in the sense that there exists one interaction with a SM fermion and a sterile neutrino. The simplified Lagrangian, focusing on interactions with the third generation only, reads where The parameter κ is fixed by the origin of the vector leptoquark and is either one or zero. Note that the leptoquark couples only to right-handed down-type quarks and right-handed sterile neutrinos. In addition, the coupling to gluons arises from the SU (3) c -charge. As usual, we assume that U is also odd under the Z 2 symmetry, otherwise the second line in Eq. (47) is not allowed. This type of interaction is possible both for Majorana and Dirac sterile neutrinos and in both cases only the right-handed components will interact.
Operator matching Straightforwardly, for processes with momentum transfer well below m U the effective four-fermion interaction reads Again, this can be rephrased as the operator O N b ≡ C 33 N d , or upon identifying Λ = m U . Hence, the effective interaction mimics a vector-mediated bb ↔ N N interaction and we can compare the scenario to a Z coupling only to b R and N R inducing the same operator.

Relic density
We implement the full model via FeynRules [55] and calculate the DM relic density after freeze-out using micrOMEGAs [30][31][32][33]. In Figure 7 we plot the m Nm U plane and fix at each point |x b | such that the observed DM relic density is obtained. We compare this to the EFT limit of the model described by the matching Eq.(50) together with Λ = m U . This means that for each contour of constant |x b | we associate a dashed EFT contour with constant |C N b | = |x b |, where the m U coordinates are given by and the fraction of the right-hand side corresponds to the curves in the center row of Figure 1. We find that as long as m U > m N , the freeze-out is nicely described by the effective operators. This is again in contrast to the Z , where already close to m Z = 2m N the EFT limit is broken. The reason is that for the t-channel mediator an on-shell annihilation has to go into a pair of mediators. The different positions of the lines for Dirac and Majorana neutrinos result from the different form of the annihilation cross sections, see Eqs. (11) and (12).
LHC and other constraints on the UV-model As before, we combine generic constraints in the EFT picture, Section 2.2, with model-specific constraints in Figure 7.
1. Perturbativity: again, we indicate in gray the region when the coupling becomes |x b | > √ 4π ≈ 3.5, where the theory becomes non-perturbative.
2. Collider production: the main production channel of our vector leptoquarks is gg → U U † , which does not depend on the coupling x b , but is entirely determined by the strong gauge coupling and the leptoquark mass m U [69]. The value of the coupling x b is only relevant in so far as it is assumed that the decay width is sufficiently small to consider the leptoquarks to be produced on-shell. The most stringent bound on m U can be derived from the leptoquark search in CMS [75]. Assuming m N = 0, the CMS analysis excludes m U 1558 GeV for κ = 0 and m U 1927 GeV for κ = 1 at 95% CL. Since a detailed exclusion contour of the m N -m U parameter space is only given for a scalar leptoquark, we project the excluded mass range of m U given in the limit of massless invisible final states onto larger DM masses. Using the scalar leptoquark plot in Ref. [75] as guiding example, this appears to hold approximately up to mass of about 500 GeV. In Figure 7, the LHC limit obtained this way is shown in red (κ = 0) and light red (κ = 1). We note that a collider study of a vector leptoquark coupling to massive neutral fermions is missing.
3. Indirect detection: as illustrated in the second row of Figure 1, current constraints from annihilation in dSphs can be used to set a lower bound on the DM mass in the case that the effective operator O N b induces the dominant contribution to DM annihilation. In Figure 7, we therefore include the lower bounds on m N including the band corresponding to the variation of the J-factors within their 68% CL intervals. Note again the different position of the indirect detection constraints for Dirac and Majorana fermions, which results from the different form of the annihilation cross sections, see Eqs. (11) and (12). In Section 2.2, we noted that the effective DMgluon coupling may be relevant for indirect detection. Considering that for a scalar sbottom-like mediator we found no notable effect in the parameter region where the observed relic density is produced (see the orange band in the center-left plot of Figure 1), we expect that also in the vector case the effect can be neglected.
4. Direct detection: in the discussion in Section 2.2 we saw that the coupling to nucleons induced by RG running of the operator O N b is not yet tested with sufficient precision to constrain the model in the EFT region, see Figure 1. While there exists no discussion of the effective DM-Higgs and DM-gluon coupling mediated by vector leptoquarks, we boldly extrapolate from the scalar case that the explicit mass suppression weakens the direct detection limits to a level where they do not affect our EFT parameter space.
Again, we can summarize Figure 7 in that there exists parameter space for a neutral fermion singlet that couples via vector leptoquarks to b-quarks, in which a consistent EFT description is possible. For Dirac DM m N > 34 GeV for κ = 0 (or m N > 52 GeV for κ = 1) and TeV-ish leptoquarks are still allowed, while for the Majorana case the sterile neutrino mass needs to exceed 77 GeV for κ = 0 and 123 GeV for κ = 1. The precision of these numbers would benefit from a dedicated evaluation of LHC constraints on vector leptoquark pair production with subsequent decays into b-quarks and massive dark fermions, as it has already done for scalar leptoquarks. Compared to a Z coupling only to b R and N R , which leads to the same effective four-fermion operator, the collider constraints are stronger in the leptoquark case.

Conclusions
Given the energy scales involved, an EFT approach to weak-scale DM is still an attractive scenario. It can, for instance, reveal tensions between freeze-out production, direct detection, and indirect detection. The problem with a pure EFT approach is that dark matter mediators can appear on their mass shell either at colliders or at some time during the thermal history of the Universe. This is why we consider effective theories at the weak scale only as representative generalizations of UV-complete models. If we force the interpretation of LHC searches into a DMEFT framework, it immediately leads to tensions with the annihilation rate or observed relic density [13].
Sterile neutrinos as WIMP dark matter offer a natural way out of this when they couple only to third-generation fermions. Starting with a pure EFT approach we have demonstrated that the relic density can be generated while all constraints from direct and indirect detection, as well as lepton flavor universality are obeyed. This holds true when we couple the heavy neutrino to tau leptons, bottom quarks, and top quarks. LHC constraints are strongest when we produce the mediator on its mass shell, so we expect them to depend on the underlying model.
We have confronted our successful νDMEFT with a set of plausible UV-complete models. We studied a (B − L) 3 gauge extension of the Standard Model as well as scalar and vector leptoquarks. In all cases, DM-mediator couplings to light fermions are absent at leading order. We found that all three models are well represented by the unifying EFT, but that LHC searches naturally distinguish between the fundamental models. In general, the LHC reach for leptoquark UV-completions is larger, because of their unavoidable coupling to gluons and the corresponding pair production process. Search limits from supersymmetric squarks usually apply with minor modifications, even though this phase space approximation is less motivated for the vector leptoquark case.

Appendices
A Details on the Z → τ τ limit in the (B − L) 3

model
In this section, we summarize how we produced the exclusion curve of Z → τ τ in Figure 5. We consider an EFT limit of the Lagrangian Eq.(19) along the following reasoning. At leading order, i.e. 0th order in and tree-level QCD, the only production channel of the Z resonance is bb → Z . Therefore one needs to consider a 5-flavor parton distribution function to evaluate the pp → Z production cross section. At the order 0 , we can ignore Z-Z mixing, and the interactions are simply determined by the (B − L) 3 charge of the respective fermions, where f = t, b, τ, ν τ , N . Given these interactions, it is straightforward to calculate the decay width of Z → f f , for f a fermion with both chiralities, in the limit m Z m f , Therefore the total decay width, considering the channels tt, bb, τ + τ − , ν τ ν τ , and N N reads Γ X Z m Z = 1 12π g 2 X (q t X ) 2 · 3 + (q b X ) 2 · 3 + (q τ X ) 2 + 1 2 (q ντ ) 2 + 1 2 (q N ) 2 = 2 9π g 2 X . (A. 3) The interaction in Eq.(A.1) along with mass terms for N and Z and the decay width in Eq.(A.3) was implemented as a Universal FeynRules Output (UFO) model file using FeynRules [55]. This model file was subsequently loaded into MadGraph [76], where we calculated the cross section for pp → τ + τ − to compare with the limits given by ATLAS [60].
B Details on Z with coupling only to b R or t R Here we discuss how we to test for constraints on a Z with effective couplings to only N and b R or t R , i.e. (B.1) In Refs. [61,62] limits on the production of a neutral scalar H 2 together with a b pair or t pair, respectively, have been presented in the form of cross section times branching into H 2 → bb or tt. We compare this product to the production cross section for pp → bbZ and pp → ttZ respectively, which we calculate using the interactions in Eq.(B.1) as before via MadGraph [76]. This has to be multiplied by the branching fraction into b or t pairs which, following Eqs.(A.2) with only right-handed quarks, is 1/4. By comparing the cross sections for m Z = m H 2 , we find in both cases that the cross section limits are not yet strong enough to probe the model for g X ≤ 3.5. We note that the tiny red area that is the limit in the bottom plot of Figure 7 has been obtained for a branching ratio in bb of 1, for a value of 1/4 it would disappear in the plot.