ournal of C osmology and A stroparticle hysics J Fermi Bubbles under Dark Matter Scrutiny Part II: Particle Physics Analysis

. The analysis of the gamma-ray photons collected by the Fermi Large Area Telescope reveals, after removal of astrophysical background, the existence of an excess towards the Galactic center. This excess peaks around few GeV, and its origin is compatible with the gamma-ray ﬂux originating from Dark Matter annihilation. In this work we take a closer look on this interpretation; we investigate which kind of Dark Matter, and which type of interactions with the Standard Model ﬁelds are able to reproduce the observed signal. The structure of the paper is twofold. In the ﬁrst part, we follow an eﬀective ﬁeld theory approach considering both fermionic and scalar Dark Matter. The computation of the relic density, the constraint imposed from the null result of direct searches, and the reliability of the eﬀective ﬁeld theory description allow us to single out only two viable dim-6 operators in the case of fermionic Dark Matter. In the second part, we analyze some concrete models. In particular, we ﬁnd that the scalar Higgs portal can provide a simple, concrete and realistic scenario able to explain the GeV excess under scrutiny. the Fermi bubbles present a fairly ﬂat (in E 2 γ d Φ /dE γ d Ω) spectrum, with a 3 7 GeV/cm A common relies on the assumption of the existence of a extended population of high-energy electrons trapped inside the emitting gamma rays via Inverse Compton Scattering (ICS) on the ambient Moreover, this hypothesis allows to establish, via synchrotron radiation in the presence of microgauss magnetic ﬁeld, an interesting spatial correlation with the WMAP haze [8] observed in the microwaves.


Introduction
The quest for the first non-gravitational evidence of Dark Matter (DM) in the Galaxy is becoming, as time goes by, one of the most intriguing and challenging tasks in astrophysics.
Paradoxically though it may seem -DM should be dark, after all -one of the most promising signals that could reveal the elusive fingerprints of DM is the analysis of the gamma-ray photons collected by the Fermi Large Area Telescope (LAT). A number of recent analysis, in particular, point towards the possible existence of a residual excess originating from the Galactic Center (GC), peaked at few GeV, and compatible with the annihilation of DM particles with a thermally averaged cross section, σv ∼ 10 −26 cm 3 s −1 [1][2][3][4][5], that is of the same order as the one required if DM is in thermal equilibrium in the early Universe.
Remarkably, a very similar excess has been found in refs. [6,7] studying the gamma-ray spectrum originating from the so-called Fermi bubbles. The name Fermi bubbles refers to a colossal pair of lobe-shaped structures, first observed in refs. [8,9], each extending tens of thousands of light-years above and below the Galactic plane, and covering more than one half of the visible sky, from the constellation of Virgo to the constellation of Grus. Invisible to the naked eye, they reveal themselves in full glory in the energy region from 300 MeV up to 300 GeV. Focusing the analysis of the corresponding energy spectrum on high latitudes -i.e. for |b| > 30 • , where b is the latitude in Galactic coordinates with b = 0 • corresponding to the In this section, we describe in detail the data used in this paper. As already mentioned in the Introduction, the analysis performed in refs. [6,7] reveals that the energy spectrum of the Fermi bubbles arises from the combination of two different components: i ) an ICS component, dominant at high latitudes, produced by a population of high-energy electrons trapped inside the bubbles, and ii ) an additional component, responsible for a bump at E γ ∼ 1 − 4 GeV, compatible with DM annihilation, and dominant at low latitudes. In ref. [7] these two component have been studied together, performing a fit of the whole Fermi bubbles energy spectrum in the energy range E γ = 0.3 − 300 GeV. In this paper, however, we are mostly interested in the analysis of the DM component. In order to simplify the analysis, as a consequence, we subtract the ICS component from the energy spectrum of the Fermi bubbles; this subtraction procedure, in fact, allows us to isolate the DM contribution in which we are interested. We follow the approach provided by ref. [6], and we present in appendix A the resulting data (collected in table 3 and table 4). As a representative example, we show in the left panel of figure 1 the energy spectrum of the Fermi bubbles after ICS subtraction in the region |b| = 10 • − 20 • where the bump at E γ ∼ 1 − 4 GeV clearly stands out. In the rest of the paper, we will refer to the spectrum of the Fermi bubbles after ICS subtraction as the "Fermi bubbles excess".
We are now in the position to fit the data describing the Fermi bubbles excess using the prompt gamma rays produced from DM annihilation. 1 Our purpose is to capture, in a completely model-independent way, the most general features of this signal in order to have a guideline for the rest of the analysis. The differential photon flux from DM annihilation is given by and r(s, θ) = (r 2 + s 2 − 2r s cos θ) 1/2 . Following refs. [6,7], we use the generalized NFW profile with an inner slope γ = 1.2, 2 where the scale radius is R s = 20 kpc, and at the location of the Sun r = 8.33 kpc the DM density ρ is normalized to 0.4 GeV/cm 3 . The values of the averaged J factor used in our analysis are listed in appendix A. We perform a chi-square analysis, and we show our results in the right panel of figure 1. For simplicity, we study all the possible two-body final 1 As noticed in refs. [6,7], the ICS photons produced by DM annihilation do not play an important role in the interpretation of the Fermi bubbles excess. 2 The standard NFW profile has an inner slope of γ = 1. σv ] obtained from the chi-square fit of the Fermi bubbles excess using the prompt gamma-ray flux produced by DM annihilation in eq. (2.1). We explore different annihilation channels, i.e. bb, cc, qq, τ + τ − . The goodness of the fit can be characterized by means of the reduced chi-square χ 2 min /n, where n ≡ N − p, N is the number of data points of the non-negative photon flux, and p = 2 is the number of the fitting parameters. The dashed red line corresponds to the value σv = 3 × 10 −26 cm 3 s −1 .
states one-by-one, i.e. assuming 100 % DM annihilation into each one of the SM channels. We find that the only annihilation channels that can fit the Fermi bubbles excess involve bb, cc, qq, τ + τ − in the final state. The photon spectrum produced by DM annihilation into final states involving light leptons, electroweak gauge bosons and the Higgs, in particular, can not reproduce the bump observed at E γ ∼ 1 − 4 GeV. The best fit values for the DM mass and the thermally averaged cross section vary from M DM ∼ 10 GeV, σv ∼ 6 × 10 −27 cm 3 s −1 (annihilation into τ + τ − ) to M DM ∼ 60, σv ∼ 2 × 10 −26 cm 3 s −1 GeV (annihilation into bb).
Let us now discuss the implication of figure 1 from the point of view of the DM relic density. The evolution of the DM density, according to the standard freeze-out scenario, is driven by the expansion of the Universe and the interactions of DM with SM particles; this picture finds a quantitative description in terms of a Boltzmann equation whose approximate solution is The value of the DM relic abundance measured by the Planck collaboration [11] is Ω DM h 2 = 0.1199 ± 0.0027. As a consequence, from the right panel of figure 1, one can naively conclude that only the bb final state has the possibility to reproduce the correct value of relic density within the confidence region obtained from the fit of the Fermi bubbles excess. This intuition is true, however, only if the thermally averaged cross section in eq. (2.1) and eq. (2.4) are the same. Considering the expansion, σv = a + bv 2 + O(v 4 ), in eq. (2.1) we have σv ≈ a, since at the present v 2 ∼ 10 −6 . In order to have the same thermally averaged cross section -4 -JCAP04(2014)020 in eq. (2.4), therefore, the s-wave must dominate over the p-wave also during the freeze-out epoch, i.e. a bv 2 , with v 2 ∼ 1/4. This request is far from obvious, and depends on the details of the interaction between DM and SM fermions. To fully understand the connection between the relic density and the fit of the Fermi bubbles excess, therefore, we are forced to abandon the model-independent perspective pursued in this section. In the next section, we shall investigate the interactions between DM and SM fermions using an effective field theory approach; we shall encounter a situation in which the p-wave is not negligible at the freeze-out epoch, thus drastically modifying the scenario depicted in the right panel of figure 1.
Before proceeding to our analysis, a caveat is mandatory. A different DM density profile other than the gNFW, used throughout this paper, results in different values for the averaged J-factor in eq. (2.1). For instance, going from the gNFW to the NFW profile, the averaged J-factor turns out to be reduced by 30 %. As a consequence, in the fit of the Fermi bubbles excess, the favored value of the thermally averaged cross section increases in order to counterbalance this effect.

Effective field theory approach
In the previous section, we have studied the Fermi bubbles excess treating the thermally averaged annihilation cross section σv as a free parameter, i.e. without detailing the interactions between DM and the SM sector. Taking one step further, in this section we explore to what extent the possibility to reproduce the signal depends on the structure of these interactions with the help of the effective field theory approach.
Pursuing this direction, we start with the following assumptions: 1. DM is a WIMP, and it is non-relativistic when it decouples from thermal equilibrium, i.e. it is cold DM; 2. DM is a singlet under the SM gauge group; 3. any new particle beyond the SM is much heavier than the DM particle; this assumption excludes resonance enhancement and/or co-annihilation processes.
In the following, we consider both fermionic (section 3.1) and complex scalar (section 3.2) DM.

Fermionic Dark Matter
In the framework delimited by our assumptions, the effective operators describing the interactions between DM and the SM particles take the following factorized form where i explicitly implies the contraction over Lorentz indices. Driven by the results obtained in section 2, we neglect heavy final states involving W ± , Z, and the Higgs boson. In section 4, we will loosen some of the assumptions for specific benchmark models. We study the following DM structures are the corresponding singlets. All in all we find in full generality 3 Vector: Pseudovector: Tensor: ) and G f TA in eq. (3.8) are CP-violating. 4 In the absence of CP violation in the DM sector, these operators are zero. We, nevertheless, include them in our analysis because of lack of knowledge in the DM sector. The mass insertion in eqs. (3.4), (3.5), (3.8) manifests the involvement of the Higgs doublet to break the chiral symmetry without violating the gauge symmetry, like the SM Yukawa couplings.
Analyzing the operators one by one, we base our study on the following four criteria.
1. We compute the gamma-ray flux originated from DM annihilation according to eq. (2.1). The analyzed operator must reproduce the Fermi bubbles excess, as described in section 2.
2. We compute the DM relic density solving numerically the Boltzmann equation, as summarized in appendix C, and we request each operator to reproduce the value observed by the Planck collaboration [11]. Note that if the value of the DM coupling corresponding to the correct relic density turns out to be larger than the one from the fit of the Fermi bubbles excess, then we will end up with too many photons from DM annihilation, thus contradicting the gamma ray data. If, on the other hand, the coupling is smaller, then the explanation of the Fermi bubbles excess can not be realized. 3 We consider only effective operators at dim-6. See ref. [13] for a recent discussion about the relevance of dim-8 operators for a TeV-scale DM particle. 4 In eq. (3.8), we have the following CP transformation property where (−1) 0 = 1 and (−1) i = −1 for i = 1, 2, 3. 3. For each operator, we impose the stringent bound on the DM-nucleon spin-independent (SI) elastic cross section from the null result of the XENON100 experiment [12]. Spindependent (SD) interactions are constrained by PICASSO [14], SIMPLE [15], ZEPLIN-II [16], and XENON100 [17] data but the corresponding bounds are in general weaker w.r.t. the SI ones. Nevertheless, we will comment on the role of SD DM searches in the context of our analysis.
4. Finally, we investigate the validity of the effective field theory description. We compare the best-fit value of the DM mass, obtained from the analysis of the Fermi bubbles excess, with a naive estimation of the cut-off scale, which is inversely related to the best-fit value of the coupling constant.
We summarize in table 1 the main properties of the effective operators in eqs. (3.4)-(3.8) associated with the DM annihilation cross section (relevant for the computation of the relic abundance and the fit of the Fermi bubbles excess), and the SI elastic cross section on nuclei relevant for DD. The explicit formulae are given in appendix B. Moreover, the color code (see the caption in table 1) is used to indicate the main reason why the corresponding operator with a certain final state fails; for instance, the tensor operator with the bb final state can simultaneously explain the Fermi bubbles excess and reproduce the correct DM relic density as shown in the right panel of figure 3, but it needs a large coupling constant (small cut-off scale) due to the weak magnetic dipole moment interaction, hence voiding the validity of the effective field theory description. According to our color code, we mark this operator in yellow.
Our main results are presented in the plane DM mass-DM coupling, [M DM , G f ]; in figure 2 and in figure 3 we show the chi-square fit of the Fermi bubbles excess together with the contour reproducing the correct DM relic density, and in figure 4 we include the latest XENON100 constraints [12].
Before discussing the results, let us notice that the scalar, pseudoscalar and pseudovector operators in eqs. From the phenomenological point of view adopted in our analysis, the difference between Majorana and Dirac DM can be described as follows. First, from eq. (2.1), there exists a factor 2 of difference in the photon flux for a given value of σv f . In addition, the cross section for Majorana DM is bigger w.r.t. the Dirac case given the same coupling constant. 5 As a result, going from Dirac to Majorana DM, the confidence region obtained from the fit of the Fermi bubbles excess will shift downward in the plane [M DM , G f ]. The similar behavior occurs to the computation of the DM relic density (see appendix C for details) so that the relative position of the confidence region obtained from the fit w.r.t. the contour reproducing the DM relic density is the same for Majorana and Dirac DM.
We summarize the results of our analysis as follows.
• Scalar operator O f S . The annihilation cross section is mass and velocity suppressed. As a consequence, due to the small DM velocity today (v ∼ 10 −3 ), this type of operators is not able to produce the Fermi bubbles excess, and henceforth it will be neglected.
• Pseudoscalar operator O f PS . We start our discussion with Dirac DM. The annihilation cross section is mass suppressed but not velocity suppressed, thus excluding immedi-  Final states involving τ + τ − , cc, bb provide a good fit to the Fermi bubbles excess, as shown in the left panel of figure 2, where we plot the 68 % and 95 % confidence regions. Only the confidence region obtained considering the bb final state leads to the correct DM relic density while cc and τ + τ − channels do not have overlap. In terms of DD searches, in the non-relativistic limit, the operator of G f PS (G f PSA ) has SI (SD) interactions [19]; the leading order term of the DM-nucleus scattering is proportional to S χ · q ( S χ · q S N · q), where S χ ( S N ) is the spin of the DM particle (nucleus) and q is the transferred momentum, typically of O(MeV) for DD experiments. Therefore, the pseudoscalar operator can not be constrained by direct search experiments because of the large momentum suppression [20]. Finally, it is worthwhile to comment on the reliability of the effective field theory approach. Considering, for example, the bb final state, the best-fit value of the pseudoscalar coupling G b PS ∼ 6 · 10 −7 GeV −3 leads to a naive estimation of the cut-off Λ ∼ (G b PS ) −1/3 ∼ 120 GeV. On the other hand, with the corresponding best-fit value of the DM mass M χ ∼ 57 GeV, we argue that the effective field theory approach based on the assumptions enumerated in section 3 might not be trustable in this case. The quantitative estimate of the validity of the effective approach can be characterized by the ratio s/Λ 2 , resulting from the expansion, where √ s 2M χ . 6 This parameter corresponds to the difference between the zeroth and first order approximation in the limit of s Λ. We present these values in the last 6 Notice that in our naive estimation of the cut-off, i.e. G f In the following, we will always adopt this assumption. However, it is straightforward to translate our results into general cases of c f PS = 1.  Table 1. Properties of scattering cross sections considering fermionic DM. For each operator we highlight the corresponding behavior concerning annihilation processes (3 td and 4 th column), and SI elastic cross sections on nuclei (5 th and 6 th column). The symbol √ (×) marks a property possessed (not possessed) by the corresponding operator. The symbol 1L indicates that the corresponding elastic cross section arises at one-loop level. We refer to the appendix B for the analytical expressions. We focus our analysis on final states of τ + τ − , bb, cc and qq (2 nd column), as suggested by the results of section 2, and we indicate in parenthesis the χ 2 min value obtained in the fit of the Fermi bubbles excess. In the last column we put the value of the ratio s/Λ 2 ; this value estimates the goodness of the effective field theory approach [see text, eq. (3.10)]. Color key. The color of each line is related to the strongest tension observed in the phenomenological analysis. Operators marked in red can not fit the Fermi bubbles excess; operators in purple are ruled out by DD experiments; operators in yellow do not give a reliable effective field theory description; operators in grey do not reproduce the observed amount of relic abundance. column of table 1. 7 We consider the effective operator approach as a good approximation whenever s/Λ 2 < 10 %. This criterium disfavors the effective theory description of annihilations into cc (43.7 %) and aforementioned bb (78.5 %). As argued before, for

JCAP04(2014)020
Majorana DM both the confidence region obtained from the fit of the Fermi bubbles excess and the contour reproducing the correct DM relic density will shift downward with their relative position fixed, thus resulting in a smaller coupling constant and a bigger cutoff scale. Therefore, the reliability of the effective approach will improve for Majorana DM. For the bb final state, however, s/Λ 2 ∼ 39.2 %, is still too large to validate the effective field theory description. To sum, the effective field theory description of DM annihilation based on a pseudoscalar operator fails to explain the Fermi bubbles excess.
• Vector operator O f V . It is well known that the annihilation cross section for the vector operator is characterized by an unsuppressed s-wave component [see eq. (B.3)], providing a natural realization of the WIMP miracle paradigm. We start with the case G f VA = 0 and present results in the right panel of figure 2. Since the DM annihilation cross section does not depend on the final state quark mass in the limit M χ m f , we have only one single orange line reproducing the correct DM relic density for the bb, cc and qq channels; the green line refers to the τ + τ − final state, which does not carry the SU(3) color charge. We find that only the bb final state can reproduce both the Fermi bubbles excess and the correct DM relic density. Additionally, for this channel, the best fit value of the coupling G GeV for the cutoff scale with s/Λ 2 ∼ 1.9 %. This is also the case for cc, qq and τ + τ − channels in the wake of the unsuppressed s-wave contribution.
We now switch to DM-nucleus scattering relevant for DD experiments. A vector operator has a non-zero SI elastic cross section on nuclei [20]; the null results of the XENON100 experiment [12], therefore, can put strong bounds on the [M χ , G f V ] parameter space. In particular, we find that the elastic cross section on nuclei mediated by the interactions of DM with the valence quarks (q=u,d) of the nucleon is ruled out by several orders of magnitude considering the XENON100 bound. On the other hand, sea quarks do not contribute at the tree level to the nuclear vector current, due to the exact cancellation between particles and antiparticles. It could lead to a superficial conclusion that the bb, cc, and τ + τ − final state are not constrained at all by the XENON100 experiment. These final states, nevertheless, can have a sizable σ SI at one-loop through electromagnetic interactions: the photon emitted from virtual fermion loop couples to the nucleus of charge Z (see figure 8). The importance of these interactions has been already emphasized in refs. [21,22] in the context of leptophillic DM (see also ref. [23]). In appendix B.4, we review the one-loop computation for channels of interest. The resulting DD bounds are presented in figure 4. Strikingly, the XENON100 null result is able to exclude a relatively large region of the parameter space even in the case of the loop-induced processes. As a result, the cc interpretation of the Fermi bubbles excess is in strong tension with the DD experimental bound (figure 4, central panel). The best fit region for the bb final state lies in the allowed region (figure 4, left panel) but well within the reach of future sensitivity [24]. Annihilation into τ + τ − final state, on the contrary, can fit the Fermi bubbles signal with light DM masses, a region where the XENON100 bound becomes weaker.
For the opposite situation, namely G f V = 0, G f VA = 0, the annihilation cross section has unsuppressed s-wave contribution, and consequently we obtain quite similar results. In particular, annihilation into bb exhibits an overlap between the Fermi bubbles excess's confidence region and the contour reproducing the correct DM relic density. The most -11 -JCAP04(2014)020 important difference comes from the DD constraint. At the tree-level, the operator (χγ µ χ)(f γ µ γ 5 f ) is velocity suppressed [19,20] and has only spin dependent interactions with nuclei. Hence it can avoid the stringent XENON100 SI bound. 8 In summary, we argue that DM annihilation into bb through the vector operator in eq. (3.6) is a forefront candidate to explain the Fermi bubbles excess in the context of the effective approach.
• Pseudovector operator O f PV . The annihilation cross section for the pseudovector operator has a s-wave contribution only in the presence of a non-zero axial-vector coupling As a result, we simply set in our analysis G f PV = 0. The s-wave is mass-suppressed due to the chirality flip, while the p-wave does not have mass-suppression. 9 This fact rules out immediately the light quark final states as explanation of the Fermi bubbles excess. Annihilations into τ + τ − , bb and cc, on the contrary, can fit the data as shown in the left panel of figure 3. As far as the computation of the relic density is concerned, only DM annihilation into τ + τ − can reproduce the observed value. Notice that this result completely overturns the scenario outlined in figure 1. The reason is that for the pseudovector operator the p-wave becomes dominant at freeze-out epoch; a naive estimation based on eq. (B.4), in fact, shows that During the freeze-out, i.e. v 2 ∼ 1/4, the p-wave dominates over the s-wave for all the final state fermions considered in our analysis. As mentioned before, this fact singles out the τ + τ − final state as the only one able to reproduce both the Fermi bubbles excess and the correct value of relic density. Moreover -analyzing the reliability of the effective field theory approach -we find that the bb and cc final states have a relatively large s/Λ 2 , while, for τ + τ − , the effective theory provides a good description. Unlike the vector operator, the pseudovector one is non-zero for Majorana DM. As argued before, for Majorana DM both the confidence region obtained from the fit of the Fermi bubbles excess and the contour reproducing the correct DM relic density will shift downward, keeping their relative position fixed. As a result, the τ + τ − channel can realize the Fermi bubbles excess while bb and cc can not reproduce the correct value 8 Considering the operator (χγ µ χ)(qγµγ 5 q) with quarks in the final state, one may be worried about the contribution to DD searches from typical anomaly triangle loop diagrams associated with two gluons. Integrating out the loop in the heavy quark limit, we find that this process is equivalent to the effective interaction where gs is the strong coupling constant andGρµ is the dual field strengthG ab = abcd G cd /2. Because of the C-and P-preserving properties of QCD, the effective operator in eq. (3.11) inherits the DD properties from (χγ µ χ)(f γµγ 5 f ), resulting in a velocity suppressed SD cross section. 9 Let us quickly recap the origin of the mass-suppressed s-wave for the pseudovector operator in eq. (3.7). Considering the DM side,χγ 0 γ 5 χ is the only component contributing to the s-wave annihilation, as can be inferred taking the non-relativistic limit [20]. In addition, considering the SM side, the combinationf γ0f is zero for a fermion-antifermion system; therefore the only possible non-zero contraction ofχγ 0 γ 5 χ involves f γ0γ 5 f . Moreover, from the CP transformation propertyψγ µ γ 5 ψ CP =⇒ (−1) µψ γ µ γ 5 ψ, it follows that both χγ 0 γ 5 χ andf γ0γ 5 f have CP = −1; this in turn implies that the s-wave annihilation must have total spin S = 0 [25] because of CP = (−1) S+1 . As a consequence, the condition S = 0 onf γ0γ 5 f forces the SM fermion-antifermion pair to carry the same helicity, implying the presence of a chirality flip (mass insertion). Following the same logic, the chirality flip is not required for the p-wave annihilation.
-12 -JCAP04(2014)020 of relic density; the only difference is that the bb and cc channels have a smaller s/Λ 2 compared to the Dirac case, 5 % for cc and 12 % for bb.
To summarize, DM (either Majorana or Dirac) annihilation into τ + τ − through the pseudovector operator in eq. (3.7) is a forefront candidate to explain the Fermi bubbles excess.
• Tensor operator O f T . The tensor structure in eq. (3.8) involves the chirality flip, and therefore leads to a mass-suppressed s-wave contribution in the annihilation cross section, as shown in eq. (B.5), excluding the light quark channels. The same behavior is shared both by the CP-preserving (G f T ) and CP-violating (G f TA ) contribution [see eq. (B.5)]. In the following, we simply set G f TA = 0; we have checked that our conclusions keep unchanged in the opposite case. The results are very similar to those of the pseudoscalar case. We plot our results in the right panel of figure 3. Final states involving τ + τ − , cc, bb provide a good fit to the Fermi bubbles excess; among them, annihilation into bb is the only channel leading to the correct amount of relic density but it suffers from the poor reliability of the effective description as a result of the weak dipole moment interactions. To sum up, the tensor operator does not provide the explanation to the Fermi bubbles excess in the context of the effective approach.
All in all, the analysis of the Fermi bubbles excess through the effective operators listed in eqs. (3.4)-(3.8) indicates some interesting scenarios. If the DM interactions with SM fermions can be described by one single effective operator, the request to reproduce the correct relic density and avoid the bound from the XENON100 experiment allows only two cases: i ) DM annihilation into bb through the vector operator in eq. (3.6), with mass M χ 52 GeV or ii ) DM annihilation (either Majorana or Dirac) into τ + τ − through the pseudovector operator in eq. (3.7), with mass M χ 10 GeV. The latter scenario relies on the interplay between the swave (dominant for the annihilation of DM today, i.e. for the computation of the gamma-ray flux) and the p-wave (compatible in magnitude with the s-wave during the freeze-out epoch).

Complex scalar Dark Matter
In this section, we study the complex scalar DM which has the following structure Similar to the fermionic DM, the effective operators in full generality can be written as Vector: Tensor:  where the antisymmetric combination ∂ [µφ ∂ ν] φ ≡ ∂ µφ ∂ ν φ − ∂ νφ ∂ µ φ preserves hermiticity. Notice, moreover, that all the terms with γ 5 are CP-violating. We analyze these operators following the same criteria adopted in the fermionic case. With the same color code employed in table 1, we summarize in table 2 the main properties of the effective operators in eqs. (3.13)-(3.16) associated with the DM annihilation cross section (relevant for the computation of the relic abundance and the fit of the Fermi bubbles excess), and the SI elastic cross section on nuclei relevant for DD. We collect the explicit formulae in appendix B. Our main results are presented in figure 5, where we show the confidence regions obtained from the fit of the Fermi bubbles excess together with the contour reproducing the correct DM relic density and the bound placed by the XENON100 experiment. We summarize the content of table 2 and figure 5 as follows. fails to provide a good and satisfactory description, introducing an error always larger than our nominal 10 % threshold of acceptance. Finally, in terms of DD searches, the scalar operator with cc and bb final states has a large SI elastic cross section. This is because, via a triangle loop, DM can interact with the gluon content of the nucleus. We briefly review the computation of the resulting cross section in appendix B.3. It turns out that the XENON100 experiment places a strong bound on the parameter space of the scalar operator. In particular, the region favored by the fit of the Fermi bubbles excess and the contour reproducing the correct relic density are ruled out.
Let us consider the opposite situation, i.e. with F s S = 0, F s SA = 0; the constraint from the XENON100 experiment goes away because the SI cross section for pseudoscalar interactions is zero. As a consequence, the annihilation into bb could provide a good fit to the Fermi bubbles excess and yield the correct relic density. 10 Notice that, however, the effective field theory approach fails to provide a good approximation due to large s/Λ 2 .
To sum, we argue that an effective field theory description of complex scalar DM annihilation based on the scalar operator in eq. (3.13) fails to explain the Fermi bubbles excess.
In section 4.1, we shall see a concrete realization beyond the effective field theory approach of eq. (3.13), in which all the tensions appearing in this analysis will be overcome.
• Vectorscalar operator O s VS . The phenomenological analysis mimics the scalar case. As far as concerns the computation of the annihilation cross section, in fact, the derivatives  on the DM side of the effective operator pick up a factor k 1 · k 2 ≈ M 2 φ at the amplitude level, where k 1,2 are the incoming DM particle four-momenta [see eq. (B.11) for the corresponding cross section]. Both the confidence region obtained from the fit of the Fermi bubbles excess and the contour reproducing the correct DM relic density are just shifted by this factor, as shown in the right panel of figure 5. We observe the same behavior also in the computation of the elastic cross section on nuclei; in this case, in fact, the kinematic of the elastic scattering implies 2k 1 · k 2 = −q 2 + 2M 2 φ ≈ 2M 2 φ , where q is the transferred momentum.
As a consequence, we reach the same conclusion as in the scalar operator: the vectorscalar operator fails to explain the Fermi bubbles excess.

Towards concrete models
The effective field theory description, discussed in section 3, is based on a set of assumptions that limit the variety of physical phenomena that can be captured. The purpose of this section is to bridge part of this gap studying two concrete models that do not fall into the territory of the effective field theory approach. In section 4.1, we study the scalar Higgs portal, 11 while in section 4.2 we propose a toy-model with fermionic DM where the mediator between the dark and the visible sector is a light leptophillic Z gauge boson.
Before proceeding, let us briefly consider some other possibilities. First, notice that a model based on the Z exchange either is excluded by the XENON100 constraint (assuming a vector-like Z-DM coupling) or can not yield the correct DM relic density (assuming a pseudovector-like Z-DM coupling). Another intriguing possibility is to realize the Fermi bubbles excess in the context of a supersymmetric model. A light Bino-like neutralino [27][28][29], annihilating into τ + τ − via the t-channel exchange of a light stau, is a potential candidate. At zero velocity, we have σv (g L g R ) 2 m 2 χ /8πm 4 τ ∼ 5.6 × 10 −27 cm 3 s −1 × (76 GeV/mτ ) 4 . As noticed in ref. [28], however, the light Bino-like neutralino can not have the right relic density without an additional entropy injection. Despite this daunting conclusion, it is still an interesting open question if supersymmetric models can accommodate the Fermi bubbles excess, and we will leave it to future study. Finally, in ref. [30] the Fermi bubbles excess has been studied in the context of the minimal hidden sector model recently proposed in ref. [31].

Higgs exchange
In this section, we add to the SM Lagrangian the following terms [32][33][34] where S is a gauge-singlet real scalar with mass, after electroweak symmetry breaking, given by m S = (µ 2 S + λv 2 /2) 1/2 . H is the SM Higgs doublet with the vacuum expectation value 11 Note added. While this work was in its final stages, ref. [26] appeared. In this paper, the Higgs portal is proposed as a possible explanation of the Fermi bubbles excess. (vev) H = v/ √ 2 = 174 GeV. The interaction term λS 2 |H| 2 is known in the literature as the Higgs portal [35]. 12 Despite its simplicity, the phenomenological implications of the renormalizable interaction in eq. (4.1) are multiple and various [35,[42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60], focusing in particular on the possibility that S can play the role of cold DM in the Universe. Imposing an extra Z 2 symmetry to ensure the stability of S, the scalar singlet can annihilate into all the SM final states due to the s-channel exchange of the Higgs boson, which could reproduce the right DM abundance (see refs. [61,62] for updated analyses). To be more quantitative, the annihilation cross section times relative velocity of the scalar singlet takes the form where Γ h ( √ s) is the off-shell decay width of the Higgs boson, including all the SM final states. We evaluate this function using the public code HDECAY [63]. 13 In the denominator of the propagator, Γ h,S is the value of the decay width of the Higgs boson at m h = 125 GeV, and it consists of two pieces: in addition to the SM contribution, Γ vis = 4.07 MeV, we have to include the invisible decay of the Higgs into two DM particles, Γ inv (h → SS), that is kinematically allowed if m S < m h /2. Therefore, we have The search for a non-zero invisible Higgs decay is currently under investigation at the LHC [64,65]. Stringent upper bounds have been placed, thus allowing to constrain, through eq. (4.5), the parameters m S and λ. Equipped with eq. (4.4), we can derive the thermally averaged cross section according to the following integral 12 It is also possible to construct a similar Lagrangian involving a gauge singlet fermionic field χ. In this case the interaction term λS 2 |H| 2 in eq. (4.1) is replaced by the following two terms [36][37][38] The first operator in eq. (4.2) leads to a velocity-suppressed annihilation cross section, and therefore it can not reproduce the Fermi bubbles excess. The CP-violating operator χγ 5 χ|H| 2 in eq. (4.2), on the contrary, is characterized by an unsuppressed s-wave. Given that in this section we are interested in minimal renormalizable realizations of the interplay between SM and the DM sector, we focus our attention on the Lagrangian in eq. (4.1). Another interesting possibility is to construct a vector Higgs portal through the interaction where Vµ is a massive gauge boson, singlet under the SM gauge group. The resulting DM phenomenology turns out to be very similar to the singlet scalar case, as discussed in refs. [39][40][41]. As a consequence, we do not dedicate a separate discussion to this realization. 13 In this way, we include both the large QCD corrections on the qq final state and the three-and four-body decays of the Higgs into W ( * ) W * , Z ( * ) Z * .

JCAP04(2014)020
where T is the Universe temperature in question and K α are the modified Bessel functions of the second kind. For the computation of the relic density, we use this definition without any approximation, solving numerically the Boltzmann equation. We briefly review the main points of this analysis in appendix C. For the computation of the photon flux, however, we use in eq. (2.1) the zero-velocity approximation σ S v = σ S v| s=4m 2 S . 14 Finally, through the Higgs couplings with quarks, the scalar singlet may also have a sizable SI elastic cross section on nuclei, σ SI , severely constrained by the XENON100 experiment. In the limit of q 2 m 2 s , where q is the momentum transferred, it has a remarkably simple form [62] where µ ≡ m N m S /(m N + m S ) is the DM-nucleon reduced mass, m N = 0.946 GeV is the nucleon mass, and f N = 0.303 is a numerical coefficient of the Higgs-nucleon coupling, which is described in the appendix B. Before checking if S as DM can explain the bubbles and avoid the experimental constraints mentioned above, we would like to point out that this simple theoretical setup should be thought as a minimal parametrization that can be used as a reference for more complicated models. For example, in the case of a complex singlet where both degrees of freedom are kinematically accessible, the relic density doubles. Even non-renormalizable derivative couplings between the singlet and the Higgs, e.g. (∂ µ |H| 2 S∂ µ S)/Λ 2 , can be constrained in a similar way. These couplings arise, for instance, in composite Higgs models where both the Higgs and the singlet are pseudo Nambu-Goldstone bosons of a broken global symmetry [67], whose breaking is characterized by the scale Λ. In this case, for an on-shell Higgs, we can identify λ → m 2 h /2Λ 2 [68]. Introducing the Higgs invisible branching ratio we fit the Fermi bubbles excess together with the results of all the Higgs searches under investigation at the LHC [69]. 15 We show our results in figure 6. The red region is excluded by XENON100, while the green line is the relic density contour; the 99 % confidence region obtained from the combined fit of the Fermi bubbles excess and Higgs data is displayed in blue. This favored region retraces the Higgs resonance peaked at m S = 62.5 GeV. The left-hand side, m S < 62.5, is cut around m S ≈ 60 GeV; for lighter DM masses, in fact, the invisible decay of the Higgs starts to be in conflict with the experimental data collected at the LHC, exceeding the allowed values. On the other hand, the region of m S > 62.5 GeV does not suffer from the LHC bound. In the region of 62.5 m S 65 GeV, the confidence region overlaps with the relic density contour while the region of m S 65 GeV is excluded by the DD bounds. 14 In the computation of the photon flux, we use the energy spectra provided by ref. [66], thus including also the effects of the three-and four-body decays h → V ( * ) V * , V = W, Z which are particularly relevant close to the Higgs resonance. 15 Performing the fit of the LHC data, we assume that the Higgs couplings with the electroweak gauge bosons and fermions are equal to their SM values. As a consequence, we have as a free parameter only, the invisible branching ratio in eq. (4.8). In this situation, ref. [69] obtains that BRinv > 22 % is excluded at 95 % C.L.. It is important to keep in mind, however, that this bound can be relaxed allowing for deviations of the one-loop Higgs couplings with photons and gluons. If so, in fact, one finds that BRinv > 50 % is excluded at 95 % C.L.. scalar Higgs portal model. We include in the fit also the LHC data constraining the invisible decay width of the Higgs, as discussed in ref. [69]. The red region is excluded at 90 % C.L. by the XENON100 experiment, while the green line is the relic density contour.
The scalar Higgs portal, in conclusion, provides a simple, concrete and realistic DM model to explain the Fermi bubbles excess without in conflict with the DD SI bounds and the LHC Higgs data. In this context, the predicted value of the DM mass lies in the small window 62.5 m S 65 GeV, a region that will be significantly reduced by the next generation of experiments (see ref. [62] for a detailed discussion).

Z -exchange
In this section, we study a toy model, where the Z boson mixes with an addition U(1) X gauge boson Z and DM is charged under U(1) X . 16 We employ the following assumptions.
1. Fermionic DM. From the effective field approach discussed above, the scalar DM annihilation via Z/Z is velocity suppressed, which cannot explain the Fermi bubbles excess.
2. Light Z . If Z is heavy, our effective field theory approach applies. In addition, the current LHC constraints [71,72] imply m Z > few TeV or m Z 100 GeV.
3. Pseudovector DM-Z interactions and leptophilic final states in order to evade the DD bound. 16 The mixing could be, for example, simply the kinetic mixing or a Stueckelberg extension of the SM [70].

JCAP04(2014)020
4. For simplicity, the couplings of Z to SM fermions are proportional to those of Z to SM fermions.
As a result, the relevant Lagrangian for Z can be written as where refers to both charged leptons and neutrinos, P L,R are the projection operators, and g R,L are the SM lepton couplings to Z. To reduce the number of free parameters, we set g χχZ = 1 and m z = 17.5 GeV. 17 Besides the Fermi bubbles excess and the DM density, we have to consider the LEP bounds. As shown in Chapter 8 of ref. [73], it can be cast into the cutoff Λ in terms of contact interactions To obtain the bound, we compute the total cross section with SM−Z interactions eq. (4.9) and compare it with the one from the contact term eq. (4.10).
In figure 7, we show a viable example of Z models, given m Z = 17.5 GeV. The green area represents the 99 % CL confidence region of the Fermi bubbles excess and the red line is the relic density contour. Similar to the Higgs exchange, the overlapping region takes place near the resonance. The LEP constraint is represented by the dashed orange line, which is well above the overlapping region between the confidence region of the fit to the Fermi bubbles excess and the relic density contour. To sum, the Z exchange with only leptonic couplings could be responsible for the Fermi bubbles excess and the correct DM density.

Conclusions
As recently revealed, the analysis of the energy spectrum of the Fermi bubbles show the existence of a component peaked at low latitude around E γ ∼ 1 − 4 GeV. In the previous work [7] -dedicated to the astrophysical analysis of this signal -we argued that its origin could be compatible with DM annihilation.
In this paper, we have analyzed this excess from the point of view of particle physics, aiming at understanding which kind of DM and which kind of interactions are able to reproduce the distinctive features of the signal without contradict the existing phenomenological constraints.
First, we have performed a model-independent analysis based on a two-dimensional fit using as free parameters the dark matter mass M DM and the thermally averaged annihilation cross section σv . We have found that the signal can be accommodated by DM annihilation into SM fermions, varying from M DM ∼ 10 GeV, σv ∼ 6 × 10 −27 cm 3 s −1 (annihilation into τ + τ − ) to M DM ∼ 60 GeV, σv ∼ 2 × 10 −26 cm 3 s −1 (annihilation into bb). LEP bound Figure 7. The Z -exchange case. The confidence region of the fit to the Fermi bubbles excess is in green and the right DM relic abundance is denoted by the red line. The Z exchange could be the underlying mechanism for both the right DM density and the Fermi bubbles excess without violating the LEP bound, the dashed orange line.
Second, we have investigated both fermionic and scalar DM using the language of effective operators in order to classify all the possible structures arising at lowest order from the interaction of DM with the SM fermions. In order to scrutinize the resulting list of operators, the Occam's razor used in our analysis relies on four basic requirements: i ) the flux of photons originating from DM annihilation must reproduce the aforementioned gamma-ray signal, ii ) the DM relic density must be in agreement with the value observed by the Planck collaboration, iii ) the SI elastic cross section on nuclei must be consistent with the stringent bound placed by the XENON100 experiment, and iv ) the best-fit values for the DM mass and coupling must validate the effective field theory description.
We have found that only two cases are able to fulfill all these four requirements: fermionic DM with mass M DM 52 GeV, annihilating into bb via a vector-type interaction, and fermionic DM with mass M DM 10 GeV, annihilating into τ + τ − via a pseudoscalartype interaction. These results can serve as guiding principles for searching for UV complete theories.
Finally, we have investigated two concrete models: the scalar Higgs portal, and a toy model of fermionic DM in which the mediator between the dark and the visible sector is a light Z gauge boson.
On the one hand, the scalar Higgs portal model turns out to be particularly suitable to accommodate all the features of the observed gamma-ray signal. This is due to the fact that when M DM 60 GeV we are close to the Higgs resonance; in this situation the dominant annihilation channel, via on-shell Higgs decay, corresponds to bb, exactly as suggested by the -21 -JCAP04(2014)020 model-independent fit performed at the beginning of our analysis. In this mass region the scalar Higgs portal can reproduce the correct amount of relic density, and avoid the bounds placed by the XENON100 and the LHC experiments.
On the other one, we have shown that a light leptophillic Z model with M DM ∼ 10 GeV and m Z ∼ 20 GeV can also provide a good phenomenological framework in which try to accommodate the observed signal without be in tension with the bound placed by the LEP experiment. In this case the dominant annihilation channel is required to be τ + τ − .
Note added in proof. After completing the paper, LUX results [74] become available, setting a WIMP-nucleon SI cross section limit at 2 × 10 −46 cm 2 . Our conclusions, however, remain mostly unchanged except for the vector operator with the bb final state in the context of fermonic DM, where part of the Fermi bubbles confidence region is in tension with the LUX results.

A Fermi bubbles energy spectrum after ICS subtraction
In this appendix, we show the data used in this paper. The data describe the energy spectrum of the Fermi bubbles after subtracting the ICS component, from an additional population of GeV-TeV electrons trapped inside the Fermi bubbles. For the subtraction procedure, we follow the same approach adopted in ref. [6]. As discussed in ref. [7], we do not use the data in the region of |b| = 1 • − 10 • , because of the large astrophysical uncertainties. We collect the energy spectrum of the Fermi bubbles after ICS subtraction, i.e. the Fermi bubbles excess analyzed throughout this paper, in table 3 and table 4. In addition, the angular average of the J factor used in section 2 for the generalized NFW profile is,J     JCAP04(2014)020 where c DM = 1 (4) for Dirac (Majorana) DM, the superscript f in G has been dropped, and we have defined the following polynomials

B.2 Annihilation cross sections: complex scalar Dark Matter
We here present for the effective operators of scalar DM the annihilation cross sections, expanded in powers of the relative velocity v where the superscript s in F has been dropped, and we have defined the following polynomials

B.3 Direct detection at the tree level
We here provide two results of the SI DM-nucleon cross section: the scalar operator for scalar DM used in figure 5, and the vector operator for fermionic DM employed to exclude the light quark final states in section 3.1.
• Scalar operator O s S . The cross section of DM scattering off quark, q, is given by • Vector operator O f V . The cross section of DM scattering off quark, q, is given by where f pu = f nd = 2, f pd = f nu = 1, and again c DM = 1 (4) for Dirac (Majorana) DM.

B.4 Direct detection at one loop: the photon exchange
We sketch in this appendix the details of our one loop computation for the elastic process describing DM scattering on a nucleus at rest with charge Z and mass m N as shown in figure 8. The differential cross-section is given by elastic cross section data, with q 2 = −2m N E nr . Throughout this work, we adopt for F(q) the Helm form factor [75]. 18 For the cross section, we finally obtain where we have defined the following form factor Notice that in logarithmical approximation thus canceling E 2 nr in the denominator of eq. (B.31) to avoid the singularity in the limit of E nr → 0.
We take the renormalization scale µ to be equal to the cut-off Λ of the effective theory, making use of the formal substitution 1/¯ + ln µ 2 = ln Λ 2 . 19

B.5 Scattering rates
In this appendix, we briefly review the basic ingredients used in the analysis of the XENON100 data. The differential scattering rate, measured in events × kg −1 × days −1 × keV −1 , is given by 18 In the present case, the Helm form factor is particularly appropriate, given that we are considering DM elastic scattering mediated by the electromagnetic interaction. Note that in general, since the DM particles are insensitive to the electromagnetic interaction, the point-like matter distribution of the proton should be used instead of the charge distributions. See ref. [76] for a recent critical discussion. 19 To be more precise, this replacement implies that the one loop dependence on the renormalization scale µ is cancelled by the same scale-dependence of the appropriate counterterms, which absorb the UV-divergences and are renormalized at the cut-off scale Λ. In the present case the counterterm is provided by the operator OM = CMχγ µ χ(∂ ν Fµν ) whose contribution to the amplitude in eq. (B. In our analysis we assume CM(Λ) = 0.

JCAP04(2014)020
where N T ≡ N Av /A is the number of nuclei in the target per unit detector mass with N Av = 6.02 · 10 26 kg −1 and A being the mass number of the target nucleus. f ( v) is the DM velocity distribution, and v min is the required minimal velocity for given nuclear recoil energy E nr . For the DM velocity distribution w.r.t. the GC, we use the Maxwellian profilẽ where N esc (z) ≡ [Erf(z) − 2ze −z 2 π 1/2 ]π 3/2 v 3 0 , v esc = 544 km s −1 , and v 0 = 220 km s −1 . In order to obtain the DM velocity distribution seen on Earth in eq. (B.38), one has to include the relative velocity of the Earth w.r.t. the GC, v obs . Therefore, we have f ( v) =f ( v + v obs ); please see ref. [77] for more details. The latest results of the XENO100 experiment [12] have been obtained analyzing 224.6 live days × 34 kg exposure. Two events have been observed in the nuclear recoil energy range E nr = 6.6−30.5 keV, consistently with the expected number of events from the background, i.e. b = 1.0 ± 0.2. In this paper we perform a chi-square analysis of the events, following the procedure outlined in refs. [78,79] (see also refs. [80,81]).

C Boltzmann equation and relic density
The expanding of the Universe and the annihilation of DM will make DM density n deviate from its equilibrium value n eq , which follows the Boltzmann equation is the Universe where a is the scale factor of the expanding Universe. For the real scalar and Majorana fermion, the symmetry factor c is 1, and for the complex scalar and Dirac fermion, c is 1 2 . We will NOT write explicitly the factor, but we include it in all the computation.
The equation can be simplified by using the variable Y = n s and x = m X T , where s is the entropy density of the Universe T is the Universe temperature, and h eff is the effective entropy degrees of freedom. The Boltzmann equation can be rewritten as The square root of the effective degrees of freedom is defined as follows, and showed in figure 9 √ g = h eff √ g eff 1 + T 3h eff dh eff dT , (C.5) where g eff is the effective energy degrees of freedom in the Universe. DM initially is in the thermal equilibrium in the Universe plasma, Y Y eq , which is the initial condition of the Boltzmann equation. In the Maxwell-Boltzmann approximation where the factor g is the internal (spin) degree freedom of the particle. At late time, after DM freezes out, its number density is much larger than its equilibrium value. Therefore, we can neglect Y eq in the eq. (C.3) and integrate the equation exactly (C.7) where the subscript f represents the value at the freeze-out time.
To obtain the DM relic density, we can solve the Boltzmann equation numerically with its initial condition eq. (C.6). Also we can solve it semi-analytically by finding the freeze-out time x f to match the initial and final value of Y eqs. (C.6), (C.7). By defining Y = (1+δ)Y eq , initially the dark matter distribution is close to Y eq and δ is small and grows slowly dδ/dt δ. In the end δ will change dramatically, since Y is close to Y f . At freeze-out time δ ∼ O(1) is required and the condition of dδ/dt δ still holds. Neglecting the derivative term of δ, the freeze-out x f can be solved by iteration of the Boltzmann equation with δ chosen to a value of O(1). Having the value of Y at present, the reduced relic density can be written as Ωh 2 2.74 × 10 8 m χ GeV Y ∞ . (C.9)