Augury of Darkness: The Low-Mass Dark Z' Portal

Dirac fermion dark matter models with heavy $Z^{\prime}$ mediators are subject to stringent constraints from spin-independent direct searches and from LHC bounds, cornering them to live near the $Z^{\prime}$ resonance. Such constraints can be relaxed, however, by turning off the vector coupling to Standard Model fermions, thus weakening direct detection bounds, or by resorting to light $Z^{\prime}$ masses, below the Z pole, to escape heavy resonance searches at the LHC. In this work we investigate both cases, as well as the applicability of our findings to Majorana dark matter. We derive collider bounds for light $Z^{\prime}$ gauge bosons using the $CL_S$ method, spin-dependent scattering limits, as well as the spin-independent scattering rate arising from the evolution of couplings between the energy scale of the mediator mass and the nuclear energy scale, and indirect detection limits. We show that such scenarios are still rather constrained by data, and that near resonance they could accommodate the gamma-ray GeV excess in the Galactic center.


Introduction
Non-baryonic dark matter (DM) accounts for about 27% of the energy budget of the universe [1]. Its particle nature is one of the most pressing puzzles at the interface of particle physics and cosmology. Several dark matter candidates have been extensively discussed and reviewed in the literature (see e.g. [2,3]); among those, Weakly Interacting Massive Particles (WIMPs) stand out for arising in several compelling particle physics models, such as supersymmetry, for naturally accounting for the DM abundance in the universe through the thermal freeze-out paradigm, and for potentially being testable with current and future experimental probes (see e.g. [4][5][6]).
The key strategies for WIMP searches are direct, indirect, and collider searches. The former consist of measuring nuclear scattering events with recoil energies on the order of the keV in underground laboratories [7][8][9][10]. WIMP signals in a direct detection experiment are directly proportional to the local dark matter density, thus the observation of a signal can be strongly tied to the presence of WIMP scattering.
Indirect detection attempts to detect the stable Standard Model particle products of dark matter annihilation, such as gamma-rays, cosmic-rays, or radiation at lower frequency in the electromagnetic spectrum [11][12][13][14][15]. The signal observed is proportional to the integrated line-of-sight dark matter density squared in the region of interest.
Finally, collider searches hinge on the fact that high-energy proton-proton collisions at the LHC can generate dark matter particles in association with other exotic particles. The associated signature would consist of missing energy in, for instance, monojet or dijet searches. Whilst not capable to unveil the astrophysical connection of the particles produced, collider studies can provide a complementary and sometimes more effective way to constrain dark matter models [16], especially with light dark matter particles.
The efficacy of each detection strategy at probing WIMPs is rather modeldependent; however, and rather interestingly, for the model we focus on this paper, there is a remarkable degree of complementarity across direct, indirect, and collider searches.
The observation of WIMP events at any of the detection strategies would be paramount to understand the laws of nature at fundamental scales, since WIMPs are expected to be embedded in UV complete models such as the minimal superymmetric standard models or minimal left-right model [17][18][19][20]. In other words, the discovery of WIMPs is tightly related to uncovering hints about underlying physics beyond the Standard Model.
In order to map the interactions between WIMPs and standard model particles which are allowed by data, simplified models have become powerful tools. In particular, simplified models which make use of vector mediators [21].
Models with a Z neutral gauge boson portal between dark and ordinary matter have attracted significant attention for a variety of reasons: they for instance represent "simplified model" version of several compelling particle models, and are constrained by data in a rather stringent way, albeit the couplings of the new boson to dark and ordinary matter are largely model-dependent [22][23][24][25][26][27][28][29][30][31][32].
Assuming the dark matter particle to be a Dirac fermion, many analysis have been done in the context of heavy mediators (M Z > 1 TeV) . The key results are that these models are plagued with restrictive spin-independent direct detection limits as well as LHC bounds on the Z mass from heavy resonance searches, limiting the allowed parameter space to the Z resonance, i.e. when the mass of the dark matter is close to half the mass of the Z .
In this work, we investigate an alternative scenario by turning off the vector coupling to Standard Model fermions as proposed in [42] to weaken direct detection bounds, and by focusing on relatively light Z masses, (M Z < 500 GeV) , to circumvent the usual heavy-resonance searches at the LHC 1 .
The present analysis markedly differs from previous analysis for a variety of reasons: (i) We focus on a very specific class of Z models, namely those where the Z possesses purely axial-vector couplings with SM fermions, and we perform a detailed dark matter phenomenology study; (ii) We show that the Z mass can be as low as 15 GeV, where the heavy resonance searches at the LHC searches are not applicable. We explicitly compute the collider limits in that region, with no rescaling, using the CL S method employing dimuon data from the LHC; (iii) We discuss the possibility of accommodating the gamma-ray excess observed in the Galactic center in the context of this class of models.
The paper is structured as follows. We introduce the model under consideration in Section 2. Section 3 is devoted to a detailed study of the invisible Z searches at LHC, whereas direct detection constraints are analyzed in Section 4. After a discussion on the Galactic center excess in Section 6, we conclude.

Model
We investigate here a U (1) X extension of the Standard Model expected to be less constrained by collider, direct and indirect detection searches. The model is based on the gauge group SU (3) c ⊗ SU (2) L ⊗ U (1) Y ⊗ U (1) X . Augmenting the SM by a new Abelian symmetry implies the existence of a new gauge boson Z , which can gain mass in different ways. To preserve gauge invariance such gauge boson will couple to SM fermion through the covariant derivativesf If q f L = q f R , i.e. the left and right-handed SM fermions transform in the same way under U (1) X (vector-like fermions), the Z will have only vectorial couplings with SM fermions, corresponding to a dark photon. Conversely, if q f L = −q f R , only axial-vector current are non-vanishing. The latter is the scenario we are interested in. The addition of a Dirac fermion dark matter field is trivial and follows the same logic. Focusing on the latter the final Lagrangian reads where χ is the dark matter candidate We remark that in order to write a Lagrangian of the form Eq. (2.2) it is necessary to assume that SM fermions be charged under the U (1) X symmetry. One should also notice that the model is clearly anomalous: due to the chirality of the SM fermions, the triangle anomalies U (1) 3 X do not cancel. Anomaly cancellation generically requires the existence of new fields. The new fields can, however, be vector-like under the SM gauge group, while being chiral under the new Abelian symmetry. With appropriate charge assignments one can construct an anomaly-free model where the Z has only axial-vector coupling to fermions. In Ref. [64], the authors have put some effort in coming up with UV complete models where the Eq. (2.2) is realized. We will thus assume that the exotic fermions needed to cancel the anomalies are sufficiently heavy so as not to spoil the dark matter phenomenology 2 . We emphasize that this assumption is crucial to the validity of our results, especially because we will be focusing on Z masses below 1 TeV.
All the numerical computations will be carried under the assumption g χv = g χa = g χ . Keeping them in the same order is arguably a natural choice. Mild departures from this assumption will change neither the relic density nor the annihilation cross section today since they are both dominated by the vectorial term. As for WIMPnucleon scattering rates, the impact is also mild. However, had we set g χv to zero, we would have been discussing a Majorana fermion, where the annihilation cross section is helicity suppressed, and the WIMP-nucleon scattering is purely spin-dependent. An overall minus sign between the couplings will induce no change to our results. Furthermore, this choice conveniently reduces the number of free parameters of the simplified model. That said, as long as one does not dramatically deviates from g χv ∼ g χa , our conclusions will readily apply.

Collider Constraints on Light Z Models
Searches for high-and low-mass dilepton resonances at the LHC have been an excellent probe of models containing new neutral vector bosons [65,66]. In the case where the new vector boson mediates the interaction between the SM and the dark sector, constraints from dijets and monojet searches for the Z are complementary in the mass versus coupling plane [36]. These are the most stringent constraints for leptophobic dark Z models. When couplings to leptons are sizable, though, dileptons searches have the potential to exclude larger portions of the models' parameter space [67][68][69] compared do dijets. This can be understood in view of the relative size of the production cross section for dijets and dileptons and their correspondent irreducible backgrounds: First, both production mechanisms are electroweak processes; second the dominant backgrounds for dijets and dileptons are the QCD jet pair production and the Drell-Yan processes, respectively. For universal fermion couplings as those assumed in this work, the relative number of flavors and color multiplicity leads to the relation (at LO) σ(pp → Z → jj)/σ(pp → Z → + − ) = 15, where denotes electrons or muons. On the other hand, at LO, for the dominant backgrounds we have σ(pp → Z → + − )/σ(pp → jj) ∼ O(10 −4 ) at the 13 TeV LHC [70], and a similar ratio should be expected at 7 and 8 TeV center-of-mass energies.

Signal simulation and branching ratios
In order to evaluate the constraints from the 7 TeV LHC data [65] below the Z pole, and above it with 8 TeV data [66], we implemented the axial Z model in FeynRules [71] to simulate our signal events. We also obtained the partial widths for the Z decays to leptons, jets, dark matter pairs and top pairs. The branching ratios and cross sections depend on four basic parameters: {M Z , M χ , g χ , g f }, the mass of the Z , the dark matter mass, the Z coupling to χ, and the (axial) Z coupling to the SM fermions, respectively.
In Fig. (1) we show the Z branching ratios as function of its mass for some benchmark points. In the upper left panel we fixed M χ = 100 GeV and g χ = g f = 0.1. We see that decays to jets dominate, followed by invisible decays, from light to heavy Z masses, while the branching ratio to leptons (electrons or muons) is of order 3%. We also observe thresholds when the vector boson is heavy enough to decay to χ and top pairs. The picture is essentially the same as either χ gets heavier or the couplings are changed but kept equal to each other, as shown at the upper right panel and the lower left panel. However, the branching ratio to dark matter reaches almost 90% when g χ g f . In this regime it is possible that a monojet search becomes as competitive as the dileptons concerning the exclusion constraints from collider data.
In the g f versus M Z plane, the branching ratio to leptons (muons or electrons) and to invisible (DM plus neutrinos) are shown in the Fig. (2). In the upper, middle, and lower rows we display the branching ratios for M χ = 10, 50, and 500 GeV, respectively. In the left(right) column we fixed g χ = 0.1(4π). The panels are split into two sub-panels: at left, the branching to leptons, and at right, to invisible.
In the weak DM-Z coupling regime (g χ = 0.1) and lighter DM masses (M χ ≤ 50 GeV), the branching ratio to electrons or muons reaches 4.5% for all g f until the top channel opens. The DM decays are low for all g f as can be seen at the right subpanels. In these scenarios, the dijet channel is the dominant one. As the DM masses increases, a heavy Z decays mainly to DM as g f gets small, reaching a 90% rate for g f ∼ 0.1. At the limit of the perturbative regime (g χ = 4π), a Z decays to DM predominantly, unless g f 0.4. The branching ratio to leptons is considerably suppressed in these scenarios, being at the 1% level for g f ∼ 1 as we see in the right column of Fig. (2).

Searches for dimuon resonances at the 7 and 8 TeV LHC
Searches for dileptons pairs with invariant masses as low as 15 GeV have been performed by the CMS collaboration [65] at the 7 TeV run with 4.5 fb −1 . Higher invariant masses up to 4.5 TeV were probed at the 8 TeV LHC by ATLAS with ∼ 20 fb −1 [66], for example, both in the dielectron as in the dimuon channel. We use the low and high mass dimuons from the CMS and ATLAS results, respectively, in order to investigate the collider constraints on the model. Signals for muon pair production were generated with MadGraph [72] with one extra QCD jet, and then interfaced with Pythia [73] for showering and hadronization simulations. Detector effects and jet clustering were taken into account with Delphes [74]. Jet matching were performed in the MLM scheme [75]. The backgrounds, as the data, were taken from the experimental studies [65,66].
The dimuons pairs were selected according to the following criteria: Low mass region  In the 15 < M < 100 GeV invariant mass region, CMS 7 TeV [65] adopted very loose criteria to select dimuon pairs: High mass region To search for high mass resonances, M > 100 GeV, with muon pairs, ATLAS 8 TeV [66] impose somewhat tighter cuts Moreover, the muons are required to be isolated. We adopted the same isolation criteria of the experimental collaborations in the Delphes settings. Be aware the slightly stronger limits are currently available from the LHC run-II with 13 TeV using 13.3 fb −1 of data for m Z > 500 GeV [76]. We estimate these limits to be stronger by a factor of 1.3 on the Z mass. Since our focus is on light Z gauge bosons, and our conclusions do not change even with the inclusion of more recent data, we simply keep this older data set.

Statistical analysis and estimated bounds
To estimate the bounds imposed on the Z masses and couplings we compared the dimuon invariant mass distributions of signal, background and data in the low and high mass regions with where d i , b i and s i represent the i-th bin count of the M distribution for data, background and signal, respectively. Our model have two free parameters: µ s for signal and µ b for the background normalization. The µ b parameter is set to the best value that fits the data for a given µ s . We employ the CL S method [77] to determine the 95% confidence limit regions on the M Z versus g f parameter space. First we calculate the related q-statistic: q(µ s ) = χ 2 (µ s ) − χ 2 (μ s ) if µ s >μ s , and 0 otherwise, whereμ s is the best fit for the signal strenght. After that we obtain the bounds by requiring The function Φ is the cumulative probability function of the standard normal distribution and q A (µ s ) is the value of the q-statistic calculated assuming d i =μ b b i , that is, when data are assumed to be represented by the best background model. Fixing the DM mass and its coupling g χ to the Z boson, we seek for the solution to Eq. (3.4) in the (M Z , g f ) plane as shown in Fig. (3).
In the upper left panel we show the m χ = 10 GeV case for three different g χ values: the lower green lines for g χ ≤ 0.1, the middle red ones for g χ = 1, and the upper black ones at the boundary of the perturbative regime g χ = 4π. The lines are discontinued at M Z = 100 GeV. The constraints for the M Z < 100 GeV were derived using the low mass region data of [65], whilst those in high mass region M Z ≥ 100 GeV with data from [66]. First, we observe that the excluded regions get larger as g χ becomes smaller once the DM cannot compete for decays with leptons and jets as can be seen at the upper row of Fig. (2). Note that the bounds saturate for g χ < 0.1.
In the low mass region, the collider constraints are as severe as in high mass region, concerning the values of g f excluded by the 7 and 8 TeV LHC, respectively, up to M Z ∼ 50 GeV. In the Z-pole region, the constraints get softened by virtue of the huge SM Z background. Also, for heavier Z bosons, the production cross sections drop fast and the top decays are turned on rendering the σ(pp → Z ) × BR(Z → µ + µ − ) very small and again escaping the collider constraints.
As χ gets heavier, the constraints become increasingly insensitive to the coupling to the Z , once the DM channel remains closed until M Z ≥ 2m χ . This can be seen in m χ = 50, 500 and 5000 GeV panels of Fig. (3). For sufficiently heavy DM or with suppressed couplings to Z , couplings between the vector mediator and SM fermions as low as ∼ 5 × 10 −3 are excluded at 95% CL for M Z ∼ 30 and 200 GeV as we observe in Fig. (3). These particular masses are a result of the trade off among the size of Z cross section, the branching ratio to leptons, and the relative distance of the Z-pole mass region.
Comparing our 95% CL limits on g f with those of Ref. [69] for the Z − Z mixing parameter , after translating their g f f Z coupling in terms of our g f , we found agreement in their order of magnitude in the small mass region. The agreement is better for larger κ which parametrizes the level of backgrounds systematics in Ref. [69]. It should be noted that the mixed Z model [69] assumes vector-axial couplings between Z and the SM fermions, but it makes a little difference concerning the collider bounds.

Dark Matter Phenomenology
In this section we compare limits from collider searches with the constraints arising from DM phenomenology. These constraints consist in the requirement of the correct DM relic density and the compatibility with limits from both Direct (DD) and Indirect (ID) DM searches. The constraints are individually briefly illustrated below.

Relic Density and Indirect Detection
The DM relic density is determined, for the range of couplings considered in our study, by the paradigm of thermal decoupling; as a consequence the experimentally favored value Ωh 2 ≈ 0.11 [1] corresponds to a suitable value of the DM thermally averged pair annihilation cross-section. The DM features two types of annihilation channels. The first is into SM fermions. The corresponding cross-section, originated by s-channel exchange of the Z , is given by: where n c = 3 (1) for annihilations to quarks (leptons), √ s is the center-of-mass energy of the collision, and Γ Z is width of the Z : An analytic expression of the thermally averaged cross-section can be obtained through the velocity expansion [67,78]: In addition, if m χ > m Z , the t-channel inducedχχ → Z Z process is kinematically allowed. The analytic expression of σ(s) is rather contrived. We will then just report the velocity expansion given by: These analytical approximations have been validated by numerically computing the thermally averaged cross-sections through the package Micromegas [79].
Few remarks are in order: (i) Notice that as long as g χv g χa the annihilation cross-section into SM fermions is s-wave dominated, with the dark matter annihilating nearly equally to all SM fermions, except for the color index, which makes the overall annihilation to be mostly into quarks; (ii) The term that goes with g 2 χv , not helicity suppressed, gives rise to a detectable indirect detection signal at Telescopes.
(iii) The term proportional to g aχ is velocity suppressed; (iv) When the annihilation into Z pairs is turned on, even the term proportional to g aχ is no longer velocity suppressed.
(v) If we had taken g χv = 0, as would occur for Majorana dark matter, the Z resonance would not have been present, since the pole (m 2 Z − 4m 2 χ ) in the numerator cancels out with the denominator.
Keeping that in mind, we have delimited the region that sets the right relic abundance as well as the indirect detection limits from the Fermi-LAT telescope from the observation of dwarf spheroidal galaxies [80] 3 . Figure 4. Results for m χ = 10 GeV and g χ = 4π, 1 and 0.1. Combined upper bounds on the model under study, in the bidimensional plane (m Z , g f ) for the assignations of the DM mass m χ and coupling g χ reported in the different panels. The black lines delimit the correct relic density parameter space. The blue, red and green regions are excluded by LHC data. The orange region represents spin-dependent PANDA-X exclusion region, whereas the dashed curve the spin-independent LUX limit, while in purple FERMI-LAT bound.

Direct Dark Matter Detection
In the case of of a Z with purely axial couplings to quarks one would expect only the spin-dependent interaction between DM and nucleons to be sizable. These are induced by the combination of the axial couplings of the Z with DM and light quarks and the corresponding cross-section is given by (we will consider only the case of scattering on neutrons since it suffers at the moment the most stringent constraints. Notice that in the case of flavor universal couplings the scattering cross sections on protons and neutrons are substantially equal.): where g ua , g da are the vector-axial couplings between the Z' and the up and down quarks respetively, which we assume to be g f according to Eq.2.2, µ χn is the WIMPnucleon reduced mass while ∆ neut q are the quark spin fractions of the neutron. We will take these to be ∆ neut u = −0.42, ∆ neut d = 0.85, ∆ neut s = −0.08 [87]. The vectorial coupling between the dark matter fermion and the Z , g χv , is completely irrelevant for the spin-dependent scattering as one can see in Eq.4.5. Although, this coupling even if negligible in the initial Lagrangian, Eq.2.2, will be non zero, at the typical energy scales of the scattering processes since they are generated through by computing the renormalization group equations (RGE) as shown in Ref. [88] so that a spin-independent cross section is actually induced with, where g uv , g dv are the vector couplings between the Z' and the up and down quarks respetively, which we are computed through RGE effects. Because of the coherent scattering produced by spin-independent WIMP-nucleon interaction, the spin-independent limits are much more restrictive than the spindependent ones, for this reason, the spin-independent scattering even if radiatively induced may provide stronger limits in certain regions of the parameter space we we will show below. For the RGE inducedg u,dv =g u,dv (µ N ), µ N ∼ 1 GeV couplings we have adopted, for simplicity, the analytical approximation provided in appendix B of [88], retaining only the dominant contribution, induced by top quark loops, present only above the EW scale, i.e. m Z m Z . For m Z < m Z the spin-dependent limits from PANDA-X are more restrictive and for this reason the spin-independent ones below the Z-pole are not shown in the figures. In the figures we have considered the most recent limits from spin-dependent limits from the PANDA-X experiment [89], spin-independent from LUX [90].
Note that had we started with a Majorana dark matter particle from the beginning, g vχ would always have vanished, and the RG running effect would have been irrelevant. In this case, only spin-dependent limits would be applicable, the dark matter relic density annihilation cross section would not significantly change, as well as the collider bounds agreeing with [91]. Altough, we have a sizable change as far as indirect dark matter detection is concerned since in the case of Majorana (or more in general only axial couplings of the DM with the Z') DM the s-wave component of the annihilation cross-section is helicity suppressed so at late times the annihilation cross-section of the DM is small.
That said, our findings are also applicable to Majorana Dark Matter, with mild quantitative changes, by simply ignoring the Fermi-LAT limits, as well as the spinindependent limits arising from the RG running and keeping the PANDA-X spindependent bounds. At the end, the model would be less constrained by data, since the spin-independent limits from LUX rule out a significant region of the parameter space. Figure 5. Results for m χ = 50 GeV and g χ = 4π, 1 and 0.1. Combined upper bounds on the model under study, in the bidimensional plane (m Z , g f ) for the assignations of the DM mass m χ and coupling g χ reported in the different panels. The black lines delimit the correct relic density parameter space. The blue, red and green regions are excluded by LHC data. The orange region represents spin-dependent PANDA-X exclusion region, whereas the dashed curve the spin-independent LUX limit, while in purple FERMI-LAT bound.

Summary of results
The results of our DM analysis are summarized in Figs. (4)(5)(6). Here we have superimposed, for the benchmarks considered in fig. (3), the collider limits from di-muon searches with the isocontours of the correct DM relic density, the limits from spindependent cross-section, as recently determined by the PANDA-X experiment [89], spin-independent cross-section, as given by LUX [90], and the most recent limits from indirect searches of DM gamma-ray signals in DSPh [80] 4 . Figure 6. Results for m χ = 500 GeV and g χ = 4π, 1 and 0.1. Combined upper bounds on the model under study, in the plane (m Z , g f ) for the a given DM mass m χ and coupling g χ , as reported in the different panels. The black lines delimit the correct relic density parameter space. The blue, red and green regions are excluded by LHC data. The orange region represents spin-dependent PANDA-X exclusion region, whereas the dashed curve the spin-independent LUX limit; finally the purple region indicates the FERMI-LAT bound.
As already indicated, despite the radiative origin, SI interaction give stronger constraints with respect to SD ones for certain Z masses. SD limits provide nevertheless a solid complement, especially at light Z masses. Direct detection limits are competitive, or even stronger that the one from LHC for g χ 1 while the latter dominate for lower values of the DM couplings. Once the FERMI exclusion limit is taken into account, the light DM benchmark, m χ = 10 GeV is completely ruled out for g f ≤ 10 −3 . Thermal DM is still in tension with ID limits for mass of 50 GeV ad exception of the pole region, m χ ∼ m Z /2, where mismatch between the annihilation cross-section at freeze-out and at present times is induced by the so called thermal broadening [94,95].
Viable thermal DM can be obtained, far from the pole region, for higher values of the mass, e.g. m χ = 500GeV, as considered in the last row of fig. (6). Notice that, with the exception of the case g χ = 0.1, there are no regions with the viable DM relic density for m χ > m Z . Indeed because of the m 2 χ /m 2 Z enhancement and of the high values of the couplings, the DM acquires a very large annihilation cross-section into Z pair as soon as this channel becomes kinematically accessible, so that its relic density is largely suppressed with respect to the experimental expectations. For this same reason, contrary to fig. (3), there are no plots relative to m χ = 5 TeV since, in this case, the DM relic density results always several order of magnitude below the correct value, for the couplings choices.
We stress that our results are also applicable to Majorana dark matter, because had we adopted a Majorana dark matter fermion the vectorial coupling g vχ would have always been zero, and the RG running effect would have been irrelevant. In this case, only spin-dependent limits would have been applicable, with mild changes to the annihilation cross section and collider bounds. As one can see from the figures, the Majorana dark matter setup has a larger region of parameter space allowed by data, if one takes a more conservative indirect detection limit from Fermi-LAT (as we discuss in the next section). In particular, if Fermi-LAT limits are weakened, for m χ = 50 GeV, g χ = 1 as displayed in Fig.5, a much larger region of the parameter yielding the right relic abundance would be allowed by data.

Galactic Center Excess
An excess in the GeV range has been observed in the Galactic center using data from the Fermi-LAT satellite [96][97][98][99][100][101][102][103][104][105][106][107]. There are several possible astrophysical explanations for, or caveats to, this excess. An attractive particle physics solution happens to be through annihilations of 30 − 60 GeV WIMPs into quarks with an annihilation cross section of 1 − 3 × 10 −26 cm 3 s/s normalized to a dark matter local density of 0.4 GeV/cm 3 , i.e. slightly below the canonical value [104]. For the light Z model discussed here, the preferred annihilation final states is mostly to quarks, and at the resonance the annihilation cross section today is in the right ballpark of e.g. the results in [108]. Thus, the model under consideration here can indeed accommodate the GeV excess. However, current constraints from the observation of Dwarf Galaxies using Fermi-LAT data place stringent limits on the annihilation cross section today into quarks [109]. Without including uncertainties in the dark matter content of dwarf galaxies, the WIMP interpretation for the GeV excess is excluded at face value. However, a recent reassessment of the J-factor from the Fermi-LAT team, taking into account systematic uncertainties in the J-factors, weakens their limits by a factor of 2-3, thus showing that there might be still a bit of room left for the WIMP-annihilation hypothesis [80]. Our model thus offers a possible dark matter interpretation for the GeV excess, as long as a conservative limit from Fermi-LAT observation is considered. 5 .

Note
Before submission of our paper we noted the work in [91] which partially overlaps with ours, but neither incorporated the spin-independent limits resulting from RG running and indirect detection limits, nor performed a detailed collider phenomenology.

Conclusions
Dirac fermion dark matter models in the context of heavy vector mediators are forced to live near the Z resonance due to the a combination of spin-independent and LHC bounds. One may switch off the Z -fermions vectorial coupling, however, as indeed occurs in some UV-complete models, and consider light Z masses to circumvent spinindependent direct detection limits and LHC bounds on heavy resonance searches.
In this work, we have demonstrated that by including the evolution of the vector coupling between the energy scale of the mediator mass and the nuclear energy scale, this coupling, which becomes non-zero, gives rise to stringent independent limits, and that by properly deriving LHC bounds on vector mediators using the CL S method, the scenario is still rather constrained by data.
Considering a variety of data, stemming from spin-independent and spin-dependent direct detection, collider, and indirect detection, we showed that only the parameter space near the Z resonance region survives, and that one could possibly accommodate the gamma-ray excess for m χ = 50 GeV. Moreover, we have discussed the applicability of our results to Majorana dark matter models. of the European Union under the Grant Agreement PITN-GA2012-316704 ("Hig-gsTools") and has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 674896.