MeV Dark Matter: Model Independent Bounds

We use the framework of dark matter effective field theories to study the complementarity of bounds for a dark matter particle with mass in the MeV range. Taking properly into account the mixing between operators induced by the renormalization group running, we impose experimental constraints coming from the CMB, BBN, LHC, LEP, direct detection experiments and meson decays. In particular, we focus on the case of a vector coupling between the dark matter and the standard model fermions, and study to which extent future experiments can hope to probe regions of parameters space which are not already ruled out by current data.


Introduction and Statement of the Problem
The nature of Dark Matter (DM) is one of the greatest puzzles of modern particle physics, as well as the nature of its interactions with the particles of the Standard Model (SM). Despite decades of experimental effort, the only interaction between DM and SM particles that has been confirmed experimentally is gravity. However, if other interactions are present, it would be desirable to have a model independent way to test how much parameter space has been actually tested and if there is room for discovery in future experiments. This can be achieved using Effective Field Theory (EFT) techniques, in which the only light degrees of freedom are the SM particles and the DM (see [1][2][3][4] for some of the early works on the subject). The advantages of this approach are clear: it is as model independent as we can get, and it relies on just a few assumptions (namely, that the New Physics (NP) mediating the DM-SM interactions is heavier than the Electroweak (EW) scale and that it respects the SM gauge symmetry). On the other hand, since all the correlations between different operators (present in any concrete model) are lost, it is usually unfeasible to perform a global analysis involving more than a few operators. Still, many informations can be obtained, and it is on this framework that we will focus.
As is well known, most of the theoretical and experimental activity over the last decades has focused on the Weakly Interactive Massive Particle (WIMP) paradigm, i.e. a DM candidate with mass in the GeV-TeV range and with typical cross sections of Electroweak (EW) size. As a matter of fact, the region currently probed by Direct Detection (DD) experiments restrict to DM masses above 5 GeV [5,6]. In addition, one should also consider bounds from indirect detection and collider experiments, and possibly see if the simple thermal freeze-out mechanism can explain the observed DM abundance. We stress that a problem may arise in considering the LHC bounds applied to the EFT operators. Given the high centre of mass energy of the LHC, for any fixed cutoff around a few TeV a part of the produced events will have an energy above the cutoff. These events fall beyond the validity of the EFT, and as such should not be used in the computation of the bounds. This has motivated the use of simplified models [7] as an useful intermediate step between the EFT and complete models. On the other hand, as shown in [8], it is possible to obtain robust collider limits if the centre of mass energy of the event is required to be below the EFT cutoff. 1 One of the most surprising results of the analysis in [8] is that, applying naive power counting to the Wilson coefficients and making a one-coupling-onescale assumption (in the sense that only one cutoff scale Λ and one coupling g * appear in all the EFT operators), only g * 2 couplings are currently probed at the LHC. The same analysis has been applied to other operators in [12,13].
Given the plethora of null results challenging the WIMP paradigm, in the last few years the interest has turned to other regions of parameter space. In particular, the MeV region has emerged as an interesting possibility, with many well motivated models (see for example the SIMP case [14,15], some models of asymmetric DM [16] and even some supersymmetric model [17,18]). The purpose of this paper is to extend the model independent EFT analysis to the case of MeV DM, highlighting the complementarity of searches and pointing out to which extent and in which cases we should expect some signal in future experiments (especially in the g * = 1 case in which no LHC limits are available). Although less explored, some constraints are already available on the MeV DM parameter space. For instance, Cosmic Microwave Background (CMB) bounds already force the s-wave annihilation cross section into SM particles to be below the thermal one for masses below 10 GeV [19][20][21][22][23], and Big Bang Nucleosynthesis (BBN) bounds may put strong constraints on the annihilation into quarks [24]. This means that, if the DM has indeed a sub-GeV mass, the thermal freezeout paradigm has to be abandoned unless the dominant annihilation channels are p-wave suppressed. Moreover, bounds coming from colliders [12,13,25], meson decays [26,27], indirect searches [28,29] and Z-physics at LEP [30] must also be considered. Finally, the MeV region can in principle be probed in the future by DD experiments measuring DM-electron scattering [31][32][33][34][35][36][37][38] and in high intensity neutrino beam facilities [39][40][41].
As can be seen, different SM particles are involved in the processes considered. As such, it looks like the only situation in which these constraints can be combined in a meaningful way to put bounds on the EFT coefficients is when the DM couples universally to all SM particles at the scale Λ where the interactions are generated. However, this is not the only possibility. As shown in [42][43][44], dimension 6 operators mix in the running between Λ and the low energy scale at which the experiments are performed. The result is that even if some operator is not present at high energy due to some unknown selection rule of the UV theory, it will be generated at low energy by the renormalization group running. Of course, how important this mixing is in imposing bounds depends crucially on the initial value of the Wilson coefficient at the scale Λ, and on the scale Λ itself. Using this information, the complementarity of bounds in DM searches for a DM mass above 10 GeV has been explored in [44] for universal coupling to quarks, to leptons and to third generation fermions and in [45] for Z models. Moreover, bounds on pure leptophilic models coming from the LHC have been analyzed in [46]. Models with the correct relic density for MeV dark matter are given, for example, in [47,48].
The paper is organized as follows. We first briefly recall the relevant operators which we will consider throughout the article. We then present current and future constraints that apply to MeV DM. Finally, we put together all the constraints taking properly into account the Renormalization Group Equations (RGE's) and show the available parameter space for some UV configuration.

DM EFT and Running
As already mentioned in the Introduction, the main hypotheses behind DM-EFT are that the only light degrees of freedom below the cutoff are the SM particles and the DM, and that at the cutoff the whole SM gauge symmetry is respected. For definitiveness, we will always take the DM to be a Dirac singlet fermion χ with mass m DM . At the scale Λ the lagrangian is given by where L SM is the SM lagrangian and the dots represent all the operators constructed out of the SM particles only (see [49] for the complete list). In writing Eq. (2.1) we are making a one-coupling-one-scale assumption (with the coupling g * = g * (Λ) defined at the scale Λ) and we restrict the sum over dimension 6 operators only. For our purposes, this is justified since operators of dimension 5, do not mix under renormalization [43] (although they generate dimension 7 operators that can mix below the EW scale [42]). At the dimension 6 level, 32 operators are present [43]. All of them can be written as the product of a DM and a SM current where J µ χ = {χγ µ χ, χγ µ γ 5 χ} and the SM currents are given by uγ µ u, dγ µ d, eγ µ e, uγ µ γ 5 u, dγ µ γ 5 d, eγ µ γ 5 e below m Z , where the first line is appropriate in the unbroken EW phase (above m Z ) while the second one is appropriate in the broken EW phase (below m Z ). Above m Z each operator appears three times, one for each generation (for simplicity, we will assume throughout the paper that the SM currents are flavor conserving), while below m Z the top quark does not appear, since it has being integrated out. Notice that since iH † ↔ Dµ H = √ g 2 +g 2 2 (h + v) 2 Z µ , this operator does not appear below m Z because both the Z and the h bosons have been integrated out.
The anomalous dimension matrices that mix the effective operators in the running above and below m Z have been computed in [42][43][44], and are independent from the form of the DM current (since the DM current is a complete SM singlet, it does not contribute to the running). In order to compute the Wilson coefficient of the various operators at a scale µ relevant for the experiments we use the public code RunDM [50]. A comment about the interpretation of the results is in order. Our exclusions are strictly valid for experiments performed with a typical energy E Λ, since for these cases the mediator can obviously be integrated out. For experiments performed with E Λ (as can be the case for LEP II or the LHC, as we will see below), the bounds should be computed keeping the mediator that generates Eq. (2.1) in the spectrum. However, such bounds are model dependent (the details of how the mediator couples to the DM are important when considering the resonant production). As shown in [8], the approach we will take gives a model independent bound even in the region in which the resonant production of the mediator is important, in the sense that we exclude a smaller region in parameter space.
In what follows, we will focus on the so called D5 operator [4], which is the product between the vector DM current and a vector SM current. In particular, we will analyze a leptophobic case, with universal coupling to quarks only, and a leptophilic case with universal coupling to leptons only, We will briefly comment on other possibilities at the end of Section 4.

Experimental bounds
In this section we present the bounds from current or past experiments that can be applied to MeV DM. We summarize in Table 1 all the experimental bounds and the operators to which they apply.

Bounds on the annihilation cross section
The DM annihilation cross section is bounded by CMB, BBN and indirect detection constraints. Self annihilation of dark matter particles may inject hadronic or electromagnetic energy in the intergalactic medium that may alter the thermal history of the Universe. Since recombination and primordial nucleosynthesis are well understood, bounds from the CMB and from BBN are in general important.
In the case of CMB, free electrons remaining after recombination can scatter off CMB photons and modify the CMB power spectrum. CMB data from WMAP and Planck set limits on the annihilation parameter P ann ≡ f (z) σv /m DM , given in terms of the thermally averaged cross section σv and the dark matter particle mass m DM . The redshift dependent efficiency function f (z) represents the amount of energy absorbed overall by the gas, and it is species dependent. The latest constraint from the Planck Collaboration is P ann < 4.1 × 10 −28 cm 3 /s/GeV at 95% C.L. [23]. The CMB bound already rules out Experiment Process Operators involved thermal s-wave annihilation cross sections for m DM 10 GeV [19][20][21][22][23]. In the future, Cosmic Variance Limited experiments have the potential to constrain P ann < 8.9 × 10 −29 cm 3 /s/GeV [19], i.e. a factor of ∼ 5 more stringent than current bounds.
In the case of MeV DM, additional care must be taken in the computation of the CMB bound because the DM pair will annihilate to mesons rather than quarks. The coupling between mesons and the DM currents has been computed in [51] in the context of Chiral Perturbation Theory (see also Appendix A for more details), and in our computation we will consider all possible decays into light mesons and light leptons. We list in Table 1 the operators involved.
In order to impose the bounds from CMB, we use the Equations in Appendix A taking the appropriate thermal average. We set f (z) = 1 for the annihilation to mesons and we take the bound on the annihilation cross section to electrons from [52]. The choice f (z) = 1 is an overestimate of the bound. It turns out, however, that even with f (z) = 1 the meson contribution to P ann is always subdominant with respect to the electron one. Therefore, in setting the limits in Sec. 4, we will consider only the annihilation to electrons.
Turning to primordial nucleosynthesis, the injection of electromagnetic or hadronic energy in the intergalactic medium can dissociate already formed nuclei or can alter the neutron/proton ratio through pion exchange. The case of sub-GeV DM has been considered in Ref. [24]. Overproduction of 3 He put bounds on the annihilation cross section into electrons, while deuterium overproduction put bounds on the χχ →bb annihilation cross section. The bound on χχ → e + e − is always weaker than the CMB bound, while for a DM mass between 4 GeV and 20 GeV the bound on σv χχ→bb is slightly stronger than the CMB one. As we are going to see, though, in the same region the bound coming from direct detection experiments is always stronger.
Concerning indirect detection searches, bounds on MeV DM coming from diffuse Xray and Gamma ray observations have been computed in [28], while more recently bounds from cosmic rays electrons and positrons have been computed in [29] (see also [53]). In the decays of diffuse X-ray and Gamma ray, model independent bounds can be put on the annihilation cross section to electrons. For m DM 30 MeV, the limit coming from INTEGRAL e COMPTEL is of order σv rel 10 −27 cm 3 /s, while for larger DM mass the bound becomes less and less stringent until it reaches the FERMI value σv rel 10 −24 cm 3 /s for m DM 1 GeV [28]. Turning to cosmic ray data, limits can be extracted from Voyager 1 and AMS-02 [29]. For masses around m DM 10 MeV, the limits are slightly more stringent than those obtained from diffuse X and Gamma ray data. Still, they are roughly an order of magnitude weaker than those obtained from CMB.

Collider constraints: LEP and the LHC
As shown in [25], mono-photon searches at LEP II can put bounds on the operators involving electrons listed in Table 1, although only the bounds on the operators (χγ µ χ)(eγ µ e) and (χγ µ γ 5 χ)(eγ µ γ 5 e) have been computed. For m DM 20 GeV, the constraint on the two operators is the same and is as strong as Λ ∼ 500 GeV for Wilson coefficients equal to 1. As explained in the introduction, to ensure the validity of the EFT in the considered events, the cut E cm = (p DM,1 + p DM,2 + p γ ) 2 < Λ should be imposed, analogous to what proposed in [8,12,13]. A complication arises, however. The monophoton data were collected with centre of mass energies scanning between 180 GeV and 209 GeV, so that it is not completely well defined which energy scale should be used in the computation of the Wilson coefficient. Since the analysis of Ref. [25] was performed supposing E cm = 200 GeV, in what follows we will simply take all the coefficients computed at a scale µ 200 GeV, and declare that the scales below this energy cannot be probed inside the validity of the EFT.
Other signatures can be better exploited at the LHC. In particular, the strongest experimental constraints come from mono-jet searches [54][55][56][57] can be used to put bounds on the operators listed in Table 1. We recast the ATLAS search [55], imposing the cut E cm < Λ, where E cm is the centre of mass energy of the process, E cm ≡ (p DM,1 + p DM,2 + p j ) 2 . The ATLAS analysis taken in consideration allows for multiple jets and the cuts require at least one jet with a p T > 120 GeV, allowing for the presence of soft and collinear jets. We implement the dimension six operators in Feynrules [58] and use MadGraph5 aMC@NLO [59] to generate events at matrix element level with the mono-jet topology. We then pass the events to PYTHIA 6 for parton showering and hadronisation [60]. In particular, we generate 200k events at parton level with 0-, 1-and 2-jets and we perform the final recast with MadAnalysis5 [61], modifying existing code [12,62]. One of the main outcomes of the cut is that for couplings g * = 1 no bound is found, while for g * = 4π we find that, for m DM 100 GeV, the region 400 GeV Λ 12 TeV is excluded. In particular, the region Λ 400 GeV is not currently probed by the LHC, not even for large couplings.
Let us conclude with some remarks on the bounds that can be extracted from Z physics. When the "Higgs portal" operators (χγ µ χ)(iH † ↔ Dµ H) and (χγ µ γ 5 χ)(iH † ↔ Dµ H) are generated at the scale Λ, they have two effects. First, they contribute with a threshold correction to the evolution of the four fermion operators (through the SM coupling between the Z boson and the SM fermions) [43]. Second, they generate a Z − χ − χ interaction that can be bound from the Z invisible width, Γ N P inv < 1.5 MeV [30]. Both effects can be used to set stringent limits on the parameter space.

Mesons decays
For MeV DM the invisible decays of mesons play an important role. In what follows, we will focus for simplicity only on the invisible decays arising at tree level. In particular, we consider the decays Υ → χχ and J/ψ → χχ [26]. 2 Since m Υ 10 GeV and m J/ψ 3 GeV, the bounds will be relevant for m DM < 5 GeV and m DM < 1.5 GeV, respectively, and the probed Λ scales will be above the Υ and the J/ψ masses. The meson masses are also the typical energy of the process, at which the relevant Wilson coefficients must be computed. The angular momentum and C/P transformation properties of the initial state determine which operators is involved in the decay process, and the different possibilities are listed in Table 1.
The 90% C.L. constraint on the branching ratio for invisible decays of Υ(1S) and J/Ψ measured by BABAR and BES are BR(Υ(1S) → invisible) < 3.0 × 10 −4 and BR(J/Ψ → invisible) < 7.2 × 10 −4 . On the other hand, the meson decays toνν via a Z boson are negligible: BR(Υ(1S) →νν) < 9.8 × 10 −6 and BR(J/Ψ →νν) < 2.77 × 10 −8 . It is therefore enough to compute the branching fraction for the bound state to decay to DM. For each process, the bounds are practically equal for all the operators involved, and can be as strong as Λ 200 GeV for the Υ → χχ decay for an UV Wilson coefficient of order unity [26,27]. Although per se the bound is not very strong (much weaker than the LEP or LHC one, for instance), it is helpful to close in a model independent way the small Λ window left open by colliders. In the future, according to [63], we can expect roughly a factor of 10 improvement in the sensitivity of BR(Υ → χχ) at Belle II, which translates in an improvement of the bound on Λ of about a factor of 2.

Direct Detection Experiments
Direct detection experiments have set stringent constraints on the dark matter nucleon scattering cross section for dark matter masses larger than ∼ 5 -10 GeV. Indeed, for spin independent scattering LUX and Xenon1T have reached a cross section limit of ∼ 10 −46 cm 2 for m DM ∼ 30 GeV at 90% C.L. [5,6]. On the other hand, the low mass region is weakly probed. Indeed for light dark matter, the fraction of initial energy transferred to the nucleus is suppressed by m DM /m N , leading to negligible recoil energy. The LUX and Xenon1T experiments probe DM scattering with nucleons only down to masses of about 5 GeV, the maximum exclusion being σ SI 10 −42 cm 2 . However, the forecasted sensitivity of SuperCDMS in the Si and Ge modes will be able to probe the MeV parameter region down to m DM 400 MeV, with sensitivity to exclude DM-nucleon cross sections down to 10 −(39÷43) cm 2 (depending on the DM mass) [64].
Although the recoil energy for a light dark matter particle scattering off nuclei is negligible, the kinetic energy involved in the process is large enough to ionize the target atom. Experiments like Xenon10 and Xenon100 could detect the ionization of a single atom [32]. The Xenon10 experiment, using only 12 days of calibration data, can weakly probe the scattering cross section of MeV dark matter on free electrons as low as ∼ 10 −38 cm 2 . Despite the weak bound, future experiments (or the analysis of data of current experiments as Xenon100 and LUX) could produce competitive limits. Moreover, different materials and process can increase the limit on DM scattering off electrons [31,[33][34][35][36][37][38]. The most promising process seems to be the DM scattering off electrons in semiconductor targets [33], which can reach a sensitivity of about σ e 10 −(43÷42) cm 2 . We refer to Table 1 for the list of the operators contributing to the DM-electron scattering and that give an unsuppressed contribution to Spin Independent (SI) direct searches. The Wilson coefficients should be computed at the scale µ 1 GeV. 3

Relic Density
Now we discuss how to obtain the correct relic density for MeV DM. As pointed out in Section 3.1, CMB bounds generically rule out a thermal annihilation s-wave cross section for DM masses m DM 10 GeV. This leaves open the possibility that either the DM is produced thermally via p-wave annihilation, or that a non-thermal production mechanism must be invoked.
In the case of the D5 operator (see the end of Sec. 2), the annihilation cross section is s-wave and the relic abundance should be produced non-thermally. Following Ref. [65], we can compute the relic abundance inside the validity of the EFT if we suppose that the reheating temperature at the end of inflation is small enough not to produce the degrees of freedom that have been integrated out, i.e. T RH < Λ. In this case, most of the DM production happens at temperatures much larger than the mass of DM or SM particles. Considering for simplicity a universal coupling to all fermions, we get [65] where T RH is the reheating temperature, T 0 = 2.7 K the present temperature and g s * the number of effective degrees of freedom in entropy. Imposing Ω DM h 2 0.12, we get that the value of Λ able to reproduce the observed relic abundance is As we are going to see in Section 4, this region of parameter space is not currently probed, and will not be probed in future experiments.

Summary of the constraints
In this section we compare all the present bounds and future sensitivities discussed in Section 3. We are interested in determining in which cases, if any, the parameter space to be probed by future experiments is already ruled out in a model independent way by current results. In the dark matter effective field theory we have two mass scales, the dark matter mass and the cut off scale Λ, and one coupling. We will present our results in a two-dimensional parameter space (m DM , Λ), fixing the effective coupling to the maximum value allowed by perturbativity, g * = 4π, and to g * = 1. For concreteness, and to avoid bounds from structure formation [66,67], we focus on the mass range 1 MeV ≤ m DM ≤ 10 GeV. 4 We consider two benchmark models for the operator O D5 introduced at the end of Section 2: universal couplings to all quarks (Sec. 4.1) and universal couplings to all leptons (Sec. 4.2). We will comment on other possibilities in Section 4.3.

Universal couplings to quarks: leptophobic case
We start considering the case in which the DM vector current couples only to the quark vector current with flavor universal couplings. At the scale Λ > m Z the effective interactions are described by while for Λ < m Z the top has to be consistently integrated out. We would therefore expect only experiments involving interactions between quarks and dark matter to contribute. However the running of the Wilson coefficient will induce also low energy couplings with leptons c ( ) V . Solving the RGE's in the leading log approximation [44,46], we get where the Heaviside function is needed when Λ is below the EW scale. This leads to the possibility to get limits on Λ from LEP and from future DM-electron scattering experiments.
In Fig. 1 we show the excluded parameter space in the (m DM , Λ) plane. In the upper (lower) panels we show the results for g * = 4π (g * = 1), while the left (right) panels show the current (future) exclusions. Let us start with g * = 4π. In the upper left panel, the large couplings at the scale Λ lead to important LHC bounds (blue region), as discussed in Sec. 3.2. Moreover, the induced coupling with electrons is also sizeable, such that the limits from LEP also apply (this is the reason why the collider bounds goes down to Λ 200 GeV). The yellow area is excluded by mesons decay. In particular, the upper limit of about 2 TeV is set by the Υ(1s) to invisible decay (hence it applies to DM masses up to 5 GeV).  In the left panels the blue, red, yellow, green and purple regions are ruled out respectively by collider (LEP and LHC), direct detection (LUX), meson decays (BaBar and BES), CMB experiments and BBN. The grey region represents the limit of validity of the EFT, Λ < 2 m DM . In the right panel the green, emerald and orange regions will be probed respectively by future CMB, DM electron scattering and direct detection (SuperCDMS) experiments, while the grey area is already excluded by current experiments. The two upper panels consider an effective coupling g * = 4π for current (left) and future (right) experiments. The lower panel shows results for an effective coupling g * = 1.
The lower limit is instead set by a combination of the bounds of the Υ(1s) and J/Ψ decays. The J/Ψ decay sets a stronger lower limit for DM masses up to 1.5 GeV, where we clearly see the threshold due to the closure of this channel. The limits from CMB are able to cover the whole range of Λ not covered by collider and meson decays, since the annihilation cross section is s-wave (see Appendix A). As expected, the direct detection bounds from LUX are relevant only for m DM 5 GeV. Concerning future experiments (upper right panel), DM-electron scattering and limits from CMB will only be able to probe a large part of the parameter space already ruled out. Interesting information will instead come from future direct detection experiments as Super-CDMS.
Turning to g * = 1, as expected the bounds are much less severe. First, there are no LHC bounds. This is due to the issue of the validity of the EFT, as discussed in Section 3.2. In addition, the g * coupling at Λ is now too small to induce a relevant coupling to electrons, in such a way that also the LEP bound is not present. The limits from mesons decay and from LUX are weaker because they just rescale with the coupling. As for the CMB limits, we see that for DM masses above 0.01 GeV, the bound is basically a rescaling of the bound in the upper panel (although, being the coupling generated through running, there are some distortions). Below this mass we suddenly lose sensitivity due to the fact that we compute the Wilson coefficient at a scale µ 1 GeV, instead of µ 2m DM . The same happens in the right panel for the region probed by DM scattering off electrons. Future direct detection experiments as SuperCDMS will set strong limits on the scale Λ for DM masses above ∼ 300 MeV. For small couplings competitive bounds may come from DM-electron experiments for the region with small DM mass and small Λ.
Let us stress that, comparing the bounds in Fig. 1 with Eq. (3.2), we see that not even the future experiments will be able to probe the region in which non-thermal relic production is effective.

Universal couplings to leptons: leptophilic case
As a second scenario, we consider a SM gauge invariant effective field theory where the dark matter current couples universally only to leptons: In this case, the running induces low energy Wilson coefficients with light quarks [44] c (u) The presence of this couplings makes possible constraints from direct detection experiments. Indeed, the contact interaction at the scale Λ do not involve light quarks and the dark matter nucleon scattering cross section comes only from radiatively induced interactions with light quarks. The same happens for meson decays. This result is visible in the top panels of Fig. 2, where we show the constraints for this scenario with g * (Λ) = 4π. The strongest limits comes from colliders (blue region), via the LEP experiment, that exclude Λ to be between ∼ 200 GeV and ∼ 6 TeV and from CMB (green), that strongly constrains the annihilation cross section to electrons. The constraints form meson decays (yellow)  are weaker, due to the fact that the couplings to light quarks arise only radiatively, and exclude Λ between ∼ 3 and ∼ 100 GeV, for dark matter masses below 5 GeV. For a dark matter heavier than 5 GeV, the strongest limits are due to LUX (red). The right panel of the first row shows the reach of future CMB (green), DM-electron scattering (emerald) and direct detection experiments (orange).
The bottom panels show the exclusions for a coupling g * = 1. With such small couplings, the running of the Wilson coefficients is not enough to set bounds from meson decays and the bounds from LUX are reduced. The absence of the meson decay limits leaves unexplored a small region between the LEP lower limit and the CMB bound. This region will be hardly covered with the next generation CMB or DM-electron scattering experiments.
As it happens in the leptophobic case, comparing the bounds in Fig. 2 with Eq. (3.2), we see that also for the leptophilic case the region in which non-thermal relic production produces the correct relic abundance is not and will not be probed by future experiments.

Other cases
Here we discuss a few other possibilities that may arise. First, we consider the situation in which the D5 operator involves a universal coupling with all the SM fermions. Since in this case all the couplings are turned on at tree level, the running will have a rather minor effect on the bounds. In fact, we have checked that the excluded regions correspond to the strongest constraints coming from the leptophobic and leptophilic case analyzed in the two previous sections: for g * = 4π, all the region below Λ 10 TeV is probed, with the bound set by the LHC limit. On the other hand, for g * = 1 the upper bound is dominated by the LEP constraint, with Λ 500 GeV excluded. Another interesting situation is given by the so called "Higgs portal". In this case, only one between the operators (χγ µ χ)(iH † ↔ Dµ H) and (χγ µ γ 5 χ)(iH † ↔ Dµ H) is turned on at the scale Λ. The coupling to fermions arise below m Z , once the Z boson is integrated out. As discussed in Section 3.2, the Higgs portal operators induce a Zχχ coupling that can be bound by the Z invisible decay width. This bound turns out to be strong, and we have checked that for g * = 4π most of the parameter space is excluded for Λ 10 TeV, while for g * = 1 the bound is relaxed to Λ 1 TeV. Moreover, this operator induces severe constraints from dark matter scattering off nuclei for m DM 5 GeV [43].

Conclusions
As more and more parameter space is ruled out by experiments without any clear signal of dark matter discovery, it is timely to explore new venues and new regions of parameter space traditionally neglected. In this paper we have analyzed the case of dark matter with mass in the MeV range, i.e. below the reach of current direct detection experiments. This region is particularly interesting since can be probed at future direct detection experiments involving the dark matter scattering off electrons. Using a model independent approach, we have added to the standard model lagrangian all the dimension 6 effective operators that can involve dark matter, and we have properly taken into account the mixing between operators induced in the renormalization group running. Our main results are summarized in Figures 1 and 2. As can be seen, large portions of parameter space are already probed in a model independent way. Although the exact value of the maximum scale Λ already excluded highly depends on the structure and the size of the UV couplings, it is clear from the plots that, under our assumptions, most of the parameter space to which future electron scattering and CMB experiments are sensitive is already ruled out. We stress that since most of the bounds involve scales below the top mass, in this case Λ should be interpreted as the mass of some mediator generating the relevant operators. Our bound applies also to this case in the limit in which we neglect effects involving the resonant production of the mediator.

Acknowledgments
We would like to thank S. Bruggisser, S. Fichet, F. Iocco, B. Kavanagah, M. Taoso and A. Urbano for valuable discussions and suggestions. This work was supported by Fundação de Amparoà Pesquisa do Estado de São Paulo (FAPESP) and Conselho Nacional de Ciência e Tecnologia (CNPq).

A Useful Formulas
We present in this Appendix some useful formula. To set the notation, we write a generic coupling between the Dark Matter and the Standard Model fermions as Following Ref. [51], the lowest order chiral perturbation theory lagrangian coupling DM to mesons is given by where as usual the mesons hermitian matrix reads while the vector spurions including the DM currents are defined as ν µ χ = diag(c V V u , c V V d , c V V s ) Λ 2 χγ µ χ + diag(c AV u , c AV d , c AV s ) Λ 2 χγ µ γ 5 χ, a µ χ = diag(c V Au , c V Ad , c V As ) Λ 2 χγ µ χ + diag(c AAu , c AAd , c AAs ) Λ 2 χγ µ γ 5 χ. (A.4) Using the previous definitions, the annihilation cross section χχ → M M into mesons M , for a vector DM current, is given by where m M is the meson mass and the relevant couplings are