# Overcoming Gamma Ray Constraints with Annihilating Dark Matter

in Milky Way Subhalos

###### Abstract

We reconsider Sommerfeld-enhanced annihilation of dark matter (DM) into leptons to explain PAMELA and Fermi electron and positron observations, in light of possible new effects from substructure. There is strong tension between getting a large enough lepton signal while respecting constraints on the fluxes of associated gamma rays; we show how DM annihilations within subhalos can get around these constraints. Specifically, if most of the observed lepton excess comes from annihilations in a nearby (within 2 kpc) subhalo along a line of sight toward the galactic center, it is possible to match both the lepton and gamma ray observations. We demonstrate that this can be achieved in a simple class of particle physics models in which the DM annihilates via a hidden leptophilic U(1) vector boson, with explicitly computed Sommerfeld enhancement factors. Gamma ray constraints on the main halo annihilations (and CMB constraints from the era of decoupling) require the annihilating component of the DM to be subdominant, of order of the total DM density.

## I Introduction

Nongravitational signals of Dark Matter (DM) have been sought after for some time now by the astrophysical and particle physics communities. At the same time results from the Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA) experiment and from the Fermi space telescope suggest a local excess positron fraction at energies above 10 GeV as well as an excess of peaking around 500 GeV. Standard cosmic ray propagation models do not account for these excesses. An attractive explanation is that a DM WIMP (weakly interacting massive particle) is present in our galaxy at large enough concentrations to self-annihilate into standard model leptons. A TeV-scale WIMP annihilating to electron-positron pairs could produce such signals. In order to be consistent with the observed relic abundance of DM, the annihilation cross-section cm s would have to be enhanced by a factor of order 100, for example by a velocity-dependent Sommerfeld enhancement.

Many authors Cirelli:2008pk (); ArkaniHamed:2008qn (); Pospelov:2008jd (); Cholis:2008qq (); Cholis:2008wq (); Meade (); Chen:2009ab (); Watson (); Kuhlen:2008aw (); Feng:2010zp () have explored this possibility, and have constrained the allowable mass versus boost factor parameter space. However these papers assume that the dominant source of indirect signals is from annihilations in the main DM halo. In a previous work Cline:2010ag () we considered the possibility adding the effects of dark matter substructure to the theoretical model and we found examples where annihilations in subhalos could provide a significant fraction of the observed lepton excesses. We showed that one could find a better overall fit to the electron-positron data from the Fermi and PAMELA experiments, and we suggested that gamma ray constraints which are now putting considerable pressure on these models could be alleviated. Our purpose in the present work was to ascertain whether this is indeed the case.

The constraints mentioned come from recent gamma ray observations of the galaxy and from Cosmic Microwave Background (CMB) measurements. As high energy electron-positron pairs are produced and diffuse throughout the galaxy, they will emit final-state radiation as well as scatter on the ambient photon field, giving rise to 1-100 GeV gamma rays that should be detectable. Given the large expected concentrations of both DM and radiation near the galactic center (GC), gamma rays from inverse Compton scattering (ICS) near the GC are particularly constraining. The Fermi Large Area Telescope (LAT) is specifically designed to detect gamma rays in this range, and its latest results have been used to rule out large regions of parameter space for annihilating WIMP models Meade (); Cirelli:2009dv (); Strumia:2010zz (); Cirelli:2010nh ().

However in this work we will show that if a sizeable proportion of the leptons from DM annihilation originate from nearby subhalos, the constraints from GC gamma rays can be relieved. Final-state (bremsstrahlung) radiation from subhalos has been examined by other authors Bovy (); Kuhlen:2009jv (); Kistler:2009xf (); Ando:2009fp (); kuhlen (); Kuhlen:2008aw (); Hutsi:2010ai (), and ref. Kuhlen:2009is () has studied the spectrum from a nearby subhalo. In this follow-up work we extend our previous findings to a prediction of the gamma ray spectrum including a full calculation of ICS radiation in the galaxy, which we compare to the full-sky data from the Fermi LAT. We include the expected contribution to the gamma ray background coming from background electrons and positrons. Using a fully-numerical approach, we find that there is less room for new contributions from the annihilation products of the DM, making the constraints on the DM models more severe. This is a serious issue even for less cuspy and cored DM profiles, that have been shown to satisfy the constraints in previous semi-analytic treatments which ignored the background gamma ray fluxes.

In our previous paper we focused on the contributions of distant subhalos to the flux of leptons at Earth. Even though these new contributions can improve the fit to the lepton data alone, here we show that they do not soften the gamma ray constraints sufficiently to be viable. Instead, we focus on the possibility that an accidentally nearby subhalo could provide the bulk of the leptonic flux. The associated gamma rays would be sufficiently hidden by strong backgrounds if this subhalo happened to lie between us and the galactic center. The effects of nearby subhalos have been previously considered by ref. Brun:2009aj (), but only allowing for purely astrophysical boost factors, due to the density of the subhalos. Here we find that velocity-dependent Sommerfeld enhancement is crucial for obtaining a positive outcome. It is precisely because of the larger boost factor available within subhalos (which have orders of magnitude smaller velocity dispersion) relative to the main halo that we are able to soften the gamma ray constraint due to the main halo near the GC, yet have a large enough lepton signal from a nearby subhalo. In addition, we must assume that the leptophilic component of the DM responsible for these processes is subdominant to the main inert (for our purposes) component, in order to sufficiently reduce the effective boost factor for annihilations in the main halo Cirelli:2010nh (). This gives rise to the interesting possibility that different kinds of DM are responsible for the cosmic ray anomalies than those which might manifest themselves in direct detection experiments.

Using a modified version of the cosmic ray propagation code GALPROP and the data from the recent Via Lactea II simulation of dark matter evolution and collapse in a Milky Way-sized galaxy, we modelled the two-dimensional axisymmetric distribution of electrons and positrons in the galaxy. These results were combined with simulated interstellar radiation field (ISRF) data in order to compute a realistic skymap of the gamma ray spectrum expected from DM annihilation in the Galaxy, which was in turn compared with a year’s worth of diffuse gamma ray observation from the Fermi LAT.

We start with a summary the cosmic ray model and results of our previous work in Section II, before discussing the relevant ICS and gamma ray physics in Section III. In Section IV we describe our methodology, and present model-independent fits to the data in several scenarios for the distribution of subhalos and the halo profiles. In particular, we show that an accidentally nearby subhalo can provide a promising loophole to the gamma ray constraints on cuspy profiles. We also predict the gamma ray flux from the subhalo, which could provide a test of the model if future measurements and understanding of backgrounds are improved. In section V we then demonstrate that the boost factors required for this scenario can be explicitly realized in a simple class of hidden sector particle physics models. We conclude with a discussion of the overall viability of this picture in section VI.

## Ii Cosmic Ray Propagation

Inside the galactic diffusion zone, particles and nuclei propagate according to the
diffusion-loss equation Strong:1998pw (), which applies to electrons and positrons as
follows:^{4}^{4}4The full transport equation also includes the effects of convection and
diffusive reacceleration, which are mainly important for the propagation of heavier species.
Here we leave these terms out for clarity, although they were included in our full
calculations with GALPROP. These are important for determining the abundance of secondary
electrons and positrons, which come from spallation and decay of various species.

(1) | |||||

denotes the particle number density per unit momentum , represents the source function, is the spatial diffusion coefficient and is the energy loss coefficient. We seek the steady-state solution of equation (1): .

Since (1) is linear, the leptons from DM annihilation travel independently in the astrophysical background. The source comes from DM annihilation which depends on the particle physics and the local density of the dark matter:

(2) |

where the prefactor is a symmetry factor for self-annihilation, is the DM energy density, cm s is the benchmark value for standard cosmology to explain the relic density of DM, and is the energy spectrum of the annihilation products. Neglecting the effect of soft photons, the spectrum can be approximated by the simple form , where is the usual Heaviside step function, and the factor arises because that the final state has two electrons or two positrons. The latter has the correct qualitative shape, and is easier to implement in GALPROP than would be a more exact spectrum. denotes the boost factor due to Sommerfeld enhancement, originating from a nonperturbative correction due to the slow () motion of the DM particles.

To simplify our analysis, we take the boost factor to be constant throughout the main halo, and tune it to provide the best possible fit to available electron and positron data. Since the Sommerfeld effect depends strongly on velocity, typical subhalos, which have a much smaller velocity dispersion, have a much higher , and we treat it as an additional free parameter. Although each subhalo has different values of , we represent the subhalo by a single average value in this first part of our analysis, where the s are treated as being uncorrelated and best fit values are sought. This is not a limitation in the case we will eventually focus upon, namely domination of the excess lepton signal by a single nearby subhalo. A further complication is that in fact has a radial dependence within each halo, because the velocity dispersion is a function of , which has been fitted by many-body simulations such as Via Lactea II kuhlen (). We will take this into account in section V.1 by averaging over the phase space of DM in the halos, in order to make contact with the results obtained in this model-independent part of our analysis.

The spatial diffusion coefficient can be parametrized as follows Simet:2009ne ():

(3) |

Two widely-used approaches exist for solving the diffusion equation in the Galaxy: semianalytic and fully numerical. We chose the latter for Galaxy-scale propagation, in part because a numerical approach allows for better control over the spatial dependence of the astrophysical input, such as energy loss due to inverse Compton scattering. GALPROP 50.1p GALPROP:web () is a publicly available software package that solves Eq. (1) with an implicit-in-time 2D or 3D Crank-Nicholson scheme. In 2D mode, it provides a () map in cylindrical coordinates of the number density of each species within the Galactic diffusion zone. To constrain the diffusion parameters, the ratio of measured secondary-to-primary species such as B/C or sub-Fe/Fe can be simulated and fit to observations. This was done to a very high degree of accuracy in Ref. Simet:2009ne (). We used results from their best fits: , and .

The full energy loss rate is due to synchrotron radiation and inverse Compton scattering:

(4) |

is the fine structure constant and is the energy density of the galactic magnetic field, for which we used the standard parametrization:

(5) |

are the energy densities of the three main components of the interstellar radiation field (ISRF): CMB radiation, thermal radiation from dust and starlight, which lie mainly in the microwave, infrared and optical regions of the electromagnetic spectrum, respectively. GALPROP uses position-dependent maps of ISRF compiled by Porter:2005qx (), rather than using a constant energy-loss coefficient computed from a local average. The latter approach (explained in section 3 of Delahaye:2008ua ()) is commonly used in the semi-analytic model. While it is indeed quite accurate when dealing with electrons from a smooth Galaxy-wide distribution of dark matter, it is an approximation that is less precise when considering the propagation into the Galaxy of electrons from DM subhalos outside of the diffusion zone. We will nonetheless make use of the semianalytic method in Section IV.2, when only local propagation will be relevant. The position dependence of the ISRF in the Galaxy is presented in Figure 1. Further details will be discussed in section III.2.

### ii.1 Via Lactea II and GALPROP

We assumed that the DM was composed of a single Dirac fermion of mass annihilating through the channel , followed by the decay , where is some dark sector gauge boson which could also be responsible for the Sommerfeld enhancement. We considered two astrophysical models for the DM distribution: a main halo-only (MH) scenario, in which only a large, spherical halo contributed annihilation products; and a subhalo (MH+SH) scenario, where the overdensities formed by DM substructure were responsible for extra annihilation of DM into electrons and positrons. In both cases, we used a spherically symmetric Einasto profile for the DM density distribution:

(6) |

is the radial coordinate from the center of the halo, is the density at , the distance at which the slope . These parameters are simply related to the radius and rotational velocity of a given subhalo as explained in Ref. kuhlen (). The shape parameter can be read off from curve-fitting the distributions from -body simulations such as vl2 (); aquarius (). It is generally taken to be around . We took kpc for the main galactic halo, with a local dark matter density GeV cm in agreement with Via Lactea II and with other recent estimates, e.g., Catena:2009mf (). It should be noted that many authors use the convention GeV cm. This leads to a factor of difference in the constraints on the annihilation cross sections, but it is of no consequence when it comes to excluding models, since constraints come from the ratio of gamma rays lepton fluxes, which both scale linearly with .

It has been argued that direct observations of rotation velocities in the Milky Way are consistent with cored DM profiles (see for example ref. Salucci:2010pz ()). Two such examples are the isothermal and Burkert Salucci:2000ps () ansatzes. The Burkert profile has been fitted to the rotation curves of galaxies other than our own, but we are not aware of references which attempt to fit the Milky Way. To allow for the alternative possibility of a cored main halo, we will therefore restrict our attention to the isothermal profile

(7) |

adopting the values kpc and GeV/cm similar to those used by ref. Cirelli:2009dv (). These values are motivated by the constraint on the observed solar density (which we take to be somewhat higher than in Cirelli:2009dv ()) and on the mass of the Galaxy with 50 kpc as determined from circular velocity measurements. However for the subhalos we will in all cases assume the Einasto form that is suggested by Via Lactea II.

Via Lactea II vl2 () was a billion-particle simulation that tracked the evolution and collapse of particles over the history of a Milky Way-sized structure. Data about the main galactic halo and the 20,047 largest subhalos that the particles (each taken to have mass 4,100 ) merged into over the course of the simulation are available to the public. While the visible galaxy is only 40 kpc across, these subhalos extend as far out as 4000 kpc from the GC. We used the Via Lactea II subhalo data as a model for substructure sourcing electrons and positrons (from DM annihilation) at the boundary of the GALPROP diffusion zone, with an overall tunable boost factor for the subhalo annihilation rate. In addition to a larger Sommerfeld enhancement from smaller velocity dispersions within each subhalo, we expect sub-substructure unresolvable from numerical simulations to give rise to further enhancement of the annihilation cross-section. Recent estimates Kuhlen:2008aw () show that such sub-subhalos alone could increase annihilation rates by as much as a factor of 10.

Electrons from an extragalactic source have a very particular density profile. While the annihilation products from the main halo follow a roughly symmetric distribution about the GC, SH electrons sourced from the diffusion zone boundary tend to form a diffuse “shell” near the edge of the diffusion zone, as illustrated in Fig. 2. Ambient radiation prevents high-energy particles from reaching the GC, trapping them near the edge of the Galaxy. The large number of subhalos combined with a large boost factor can allow some particles to make their way to earth, albeit with a fraction of their initial energy.

We compared the best-fit combination of DM mass and boost factor for the MH scenario with the best fits for the MH+SH scenario in Cline:2010ag (). The results are summarized in table 1: a much better fit could be obtained by including subhalos and a dark matter particle with TeV, rather than the standard MH-only TeV. Of course, the fits are further improved by allowing the normalizations of the background electrons and positrons to be additional free parameters, denoted as the “freely varying background,” as opposed to the standard backgrounds resulting from GALPROP simulations which include the effects of heavier nuclear species. Assuming this extra freedom has been advocated or used by numerous authors Meade (); Cirelli:2009dv (); Cholis:2008wq (); Papucci:2009gd (). In table 1 we also show the fit we obtain in the present analysis for the main-halo-only case with an isothermal profile and fixed background. It is significantly worse than the corresponding one for an Einasto profile.

### ii.2 Annihilation channels

While we have mostly focused on the 4e final state, there is no reason for other, heavier particles not to be produced if the mass of the intermediate gauge boson is large enough. Since the amount of Sommerfeld enhancement ultimately depends on this mass, it is important to include the decays to muons and pions. The possible final states are all the four-particle combinations of , and . The muon and pion spectra are given by Ref. Cholis:2008vb (), whose authors were kind enough to provide us with the appropriate GALPROP implementation.

The branching ratios are given by , where the are given by

(8) |

In each , the square root factor comes from the phase space, while the rest is from the squared matrix element for the decay. Below threshold, is defined to be zero. For a gauge boson with a mass GeV, we find and . In this case the electrons produced from the final decay of the ’s and ’s peak at a lower energy, thus requiring a slightly higher mass of TeV in order to fit the Fermi and PAMELA data. This is much smaller than the well-known 2.2 TeV best fit in the pure-muon final state Meade (); Papucci:2009gd (); Cholis:2008wq () because of the large fraction of gauge bosons still decaying directly to high-energy electrons. These results are also shown in Table 1.

Freely-varying background (Einasto) | ||||||

(TeV) | ||||||

MH () | 0.85 | 15.5 | 18.7 | 34.3 | 90.3 | |

MH+SH | 1.2 | 2.3 | 14.2 | 16.5 | 92.8 | 3774 |

Fixed GALPROP background (Einasto) | ||||||

MH () | 1.0 | 8.2 | 144 | 152 | 110 | |

MH+SH | 2.2 | 2.1 | 175 | 177 | 146 | 1946 |

MH () | 1.2 | 3.8 | 109 | 112 | 118 | |

Isothermal profile (fixed background) | ||||||

MH () | 1.0 | 9.1 | 186 | 195 | 113 | |

MH () | 1.2 | 3.0 | 151 | 154 | 119 |

## Iii Gamma Ray Computation from Inverse Compton Scattering and Bremsstrahlung

### iii.1 “Prompt” gamma ray emission (bremsstrahlung)

Prompt gamma ray emission appears in the final stage of DM annihilation, softening the lepton spectrum. The flux can be divided into main halo and subhalo parts:

(9) |

The astrophysical and particle physics dependences of each flux can be factorized as

(10) |

and

(11) |

In each case, the factor depends only upon astrophysical inputs. The main halo factor is defined as a line of sight (l.o.s.) integral of flux at each pixel:

(12) |

In the case of flux originating from many distant subhalos, we may treat each one as a point source of radiation. In this case, the diffuse flux per solid angle requires a sum over each contributing source with density and distance within the observed solid angular region :

(13) |

This clearly depends not only on the density profiles, but also on the distribution of subhalos in the Galaxy. We will not present the results of the disant subhalo calculation of final-state radiation here, since it has been thoroughly explored by other authors in similar contexts. We direct the interested reader to references Kuhlen:2008aw (); kuhlen (); Anderson:2010df ().

Finally, if a particular subhalo is close enough to subtend an angle larger than the detector’s pixel size, it can no longer be treated as a point source: eq. (12) must be used, including the angular dependence of the projected density profile of the given subhalo, . We will return to this case in Section IV.2.

In the case of a two-lepton final state Birkedal:2005ep ():

(15) |

where and is the standard Mandelstam variable. We are interested in the case of TeV dark matter annihilating to a four-lepton final state, with a GeV leptophilic gauge boson as the messenger. The annihilation is dominated by , where the ’s are on shell. The cross section can be obtained by first computing in the rest frame of the using the decay and then boosting to the lab frame, in which the slowly moving DM particles are approximately at rest. This can easily be done numerically. We present the resulting spectrum in fig. 3. Since we will not make use of the final-state bremsstrahlung for other annihilation channels (4 or 4) we will not discuss their spectra.

### iii.2 Inverse Compton Scattering

Charged particles travelling through the interstellar medium scatter off ambient photons of the interstellar radiation field (ISRF), which is composed of microwave ( eV) radiation from the cosmic microwave background (CMB), infrared ( eV) radiation from dust, and optical ( eV) photons from starlight. Along with the galactic magnetic fields, this is the main source of energy loss for electrons diffusing within the Galaxy. We will show that ISRF photons that have scattered with TeV-scale electrons have spectra that peak at several hundred GeV, which should fall squarely within the measurement window of diffuse gamma rays by the Fermi Large Area Telescope (LAT).

Once integrated over scattering angles, the well-known Klein-Nishina formula for the Compton scattering process e e can be integrated along the line of sight to give the total flux of scattered photons per solid angle arriving on a detector Blumenthal:1970gc (); Meade ():

(16) |

represents the line-of-sight integral from the observer’s position to infinity (practically speaking, to the edge of the diffusion zone). We have used the definitions:

(17) |

and

(18) |

We numerically integrated eq. (16) along the line of sight, as well as over the incoming particle energies. All the quantities in the integrand are known: we used the two-dimensional (,) distribution of electrons and positrons from DM annihilations produced with GALPROP, as discussed in Section II. For the ISRF, we used a realistic two-dimensional photon energy density distribution from Porter:2005qx (), which is publicly available on the GALPROP website. Both distributions assumed cylindrical symmetry around the Galactic axis. For each galactic latitude-longitude pair, the line of sight integration was performed in a three-dimensional sky from the Sun’s position to the edge of the diffusion zone which was taken to extend to a radius kpc and to a height kpc above and below the galactic plane. A trapezoidal integration step size of 0.1 kpc was found to be numerically converged. The values of and at each step were found in the heliocentric coordinate system by using a bilinear interpolation scheme. On top of the DM annihilation products, we used the densities of primary and secondary electrons as well as secondary positrons to compute the ICS contribution of the background lepton field. This had the effect of further constraining the gamma ray background.

We performed the integration once per grid point on an equally-spaced latitude-longitude grid of the quarter-sky in the ranges , . This was sufficient to reconstruct the entire sky, given the symmetry of the data input.

### iii.3 Fermi all-sky diffuse gamma ray measurements

The Fermi Large Area Telescope (LAT) is a high-sensitivity gamma ray instrument capable of detecting photons in the MeV to 300 GeV range. It has an effective detector area of cm, a 2.4 sr field of view and can resolve the angle of an incident photon to 0.15 at energies above 10 GeV. Data from the first year of observation are publicly available from the Fermi collaboration.

We used the all-sky diffuse photon file from the Fermi weekly LAT event data webpage fermisky (). This covered observations from mission elapsed time (MET) 239557417 to MET 272868753 (seconds), corresponding to 55 weeks of observation between August 8 2008 and August 25 2009. We processed the photon data with the Fermi LAT science tool software, available from the Fermi Science Support Center (FSSC) website. We first removed all events with a zenith angle greater than 105 to eliminate Earth albedo. The data were further trimmed to keep only the photons measured during “good” time intervals. We then created an exposure cube from the spacecraft data for the corresponding period, to account for effective instrument exposure. The data were separated into 0.25 latitude and longitude bins spanning the entire sky, and into 16 logarithmically separated bins from 100 MeV to 200 GeV. Uncertainties were assigned according to Ref. Porter:2009sg (). We compared our results to the August-December 2008 spectrum presented by the Fermi collaboration Porter:2009sg (). The half-year data agreed exactly, while adding the extra 8 months to the full 55-week dataset changed the picture only very slightly. We rebinned the data into a 4040 grid, in correspondance with the ICS computation.

Before proceeding to the results of our numerical analyses, we should note that many factors contribute to the theoretical uncertainty. While we were able to reproduce the results of Simet et al. Simet:2009ne () quite closely, there are substantial discrepancies between the results of GALPROP and other methods of solving the transport equation. This lack of agreement is further discussed in Cline:2010ag (). There is an additional uncertainty in the injection spectrum of primary electrons, which serve, along with secondary electrons and positrons from spallation, as the astrophysical background to our results.

## Iv Empirical fits

As expected, we found that allowing subhalos to contribute to the overall flux of DM annihilation products reduced the flux of expected gamma rays from the galactic center, while increasing fluxes at higher galactic latitudes. The most stringent constraints were from the low-longitude regions just above and just below the galactic plane, where astrophysical sources of gamma rays are less prominent, but the DM distribution is still quite dense. Specifically, we used the lower right-hand region (, in Galactic coordinates) which was found to be the most constraining, in agreement with Ref. Papucci:2009gd ().

After including the ICS from background electrons and positrons, we found that the boost factor of a main halo 1 TeV DM annihilation process cannot violate the bound if the signal is to remain below the top Fermi LAT error bars. If we extend the constraint to Fermi , this condition is only slightly relaxed to . In the case of a 2.2 TeV DM candidate, these bounds become and at and , respectively. While this agrees qualitatively with other works Meade (); Papucci:2009gd (), we attribute our more stringent upper bounds mainly to our higher , as discussed in section II.1, to our inclusion of the ICS contribution from background electrons and positrons,but mainly to the different method used to solve the diffusion equation (1).

Using the best fit scenario of Ref. Cline:2010ag (), the reduction of flux was however not enough to overcome the constraints from the Fermi observations. This is illustrated in figure 4, which shows that the MH+SH scenario still violates constraints from the data by as much as . On its own, the predicted flux exceeded the data at energies above 100 GeV by at least , while we expect that additional constraints from decays should also be large in this energy range Porter:2008ve () and push predictions from this model even farther outside of the observationally allowed region. Allowing the background to freely vary (top section of Table 1) made no appreciable difference with respect to gamma rays, and was not enough to satisfy the observational constraints.

Figure 5 illustrates how the ICS gamma ray flux is increased at higher galactic latitudes when subhalos are included. It should however be emphasized that the predicted fluxes in this region of the sky are still well below the level of Fermi observations.

### iv.1 Less cuspy dark matter profiles

In section II.1 we mentioned the motivations for considering less cuspy DM profiles. Many previous works studying the ICS constraints have compared the effects of cored versus cuspy DM profiles, noting that the constraints are weaker for cored profiles. To better quantify exactly how much cuspiness can be tolerated, it is interesting to vary the parameters of the Einasto profile that control this Cholis:2009gv (); Cirelli:2010nh (). In particular, larger values of and correspond to less concentrated halos. We ran simulations of the lepton distribution and gamma ray fluxes with slightly different parameters for equation (6) while keeping the local density constant at GeV cm. This is illustrated in fig. 6. Flatter profiles with or , kpc reduce the gamma ray fluxes somewhat, but not enough to bring the predicted flux to within the observations in the offending energy bins between 10 and 100 GeV. The same is true for the isothermal profile, whose corresponding results are shown in fig. 7. For both cases, the problem arises because the predicted background gamma flux is not far below the observed flux in the most constraining bins. This leaves very little room for the additional contribution from the DM decay products ICS signal.

Increasing the intermediate gauge boson mass to 1 GeV, and thus allowing a decay to muons and pions according to the branching ratios described in Section II.2 does not alleviate the problem. Indeed, the 1 (2) bounds become BF () for an Einasto profile, and BF () in the isothermal case. These fall well short of the required BF = 118 to explain the Fermi and PAMELA excesses, as long as the DM mass is increased to TeV. These results are summarized in the bottom of Table 2. The reason ICS constraints are stronger when muons are included is due to the nature of the data. Indeed, the peak of the ICS spectrum lines up with the most constraining data point when TeV. This provides a stronger than expected constraint, relative to the 4e final state at TeV.

Subhalo | (kpc) | (pc) | (km/s) | ||

1 | 0.01 | 69 | 4.74 | 33.9 | 2.9 |

2 | 0.1 | 3.46 | 4.34 | 95.5 | 6.7 |

3 | 3.2 | 0.04 | 3.76 | 178 | 22 |

4 | 0.9 | 1.27 | 2.35 | 165 | 36 |

5 | 1.1 | 2.0 | 1.70 | 170 | 55 |

Main halo, 4e channel | |||||

Einasto | 25 | 0.048 | |||

Isothermal | 3.2 | 2.32 | |||

Main halo, 4e + 4 + 4 channel | |||||

Einasto | 25 | 0.048 | |||

Isothermal | 3.2 | 2.32 |

### iv.2 Close subhalo

The above analyses implicitly assume that no single subhalo dominates the lepton signal. But if a subhalo happens to be very close (within a kpc) to the solar system, the picture changes significantly, since the electrons and positrons from the close subhalo can dominate the observed flux, and its gamma ray emissions can come from a sizable solid angle in the sky. We treat this case separately from the previous subhalo scenario, since a larger DM mass is no longer required to produce the observed lepton signal; rather, the small amount of ICS energy loss during propagation from a local subhalo means that a 1 TeV-scale DM particle appropriately conforms to the Fermi measurements. We concentrate on the 4e final state channel, although previous results allow this to be generalized. The solution depends linearly on the spectrum , so that the boost factor required to explain the observed lepton excess should scale in the same way that it does in the main halo scenario: , as read from Table 1.

Since GALPROP is not easily adapted in its 2D mode to include the effects of a highly localized additional source term, we adopt a semi-analytic approach to solve the diffusion equation (1) for leptons produced in the nearby subhalo. Given that the leptons and gamma rays in this scenario would be from a local origin, the spatial dependence of the interstellar radiation and magnetic fields becomes much less important. We used the method described in ref. Baltz:1998xv (), with the same diffusion parameters as presented in section II (of the present work), but with an energy-loss coefficient parametrized by

(19) |

with s GeV characterizing the local energy loss rate.

We sampled subhalos from the Via Lactea II simulation to identify examples that could allow for simultaneously fitting the PAMELA/Fermi lepton fluxes and the Fermi gamma ray fluxes. Four such examples are labeled as SH1-SH4 in table 2, and a fifth (SH5) is one that we have “engineered” by choosing parameters that are close to those of SH4, but with a higher density and hence higher circular velocity, dynamically related to each other by eq. (13) of kuhlen (),

(20) |

with . Due to the higher density, SH5 requires a lower boost factor to produce the observed lepton signal, and so it represents a kind of best-case scenario. The distribution of Via Lactea II subhalos in the space of is shown as a scatter plot in fig. 8, and the five subhalos of interest are highlighted on this plot. They are atypical in the sense of needing a higher-than-average central density. A further caveat is that such a large is unlikely at small distances from the GC due to tidal disruption. Indeed, subhalos within the visible galaxy in the Via Lactea IIsimulation were of the order kpc, falling below the kpc compatible with the most plausible particle physics scenario discussed in Section V.

Each subhalo was situated along an optimal axis, namely that connecting the earth to the GC. Such an accidental alignment makes it easier to “hide” the gamma rays originating from the subhalo since they are coming primarily from the same direction as the GC, where the background emissions are strongest. This is also the reason that the most stringent ICS constraints on the main halo arise from the regions of galactic longitude instead of the most central region. However in this case we find that the biggest contribution to the emission is from final-state bremsstrahlung rather than ICS. The latter is found to produce gamma ray fluxes that are 3 orders of magnitude smaller than observed. This is consistent with the fact that the main source of ICS is IR radiation and starlight, which is concentrated far from the vicinity of the solar system.

Results were then compared to the Fermi lepton and gamma ray data in order to establish constraints. The strictest gamma constraints were at the largest energy data point from the Fermi LAT analysis of GeV, because of the shape of the FSR spectrum, which rises steadily until 1 TeV. We used a slightly different region of the sky than in our previous ICS analysis, , , because there were not enough good data points in this energy bin at lower longitudes to constrain the data. We compared the lepton prediction to the Fermi measurements at 559 GeV, where the observed spectrum is at a maximum deviation from a power law. In both cases we included the additional constraints from astrophysical backgrounds computed by GALPROP and by our ICS routine.

Results are shown in fig. 9. If the single subhalo is allowed to saturate the observed lepton signal, fig. 9 gives clear bounds (summarized in Table 2) on the proximity of each subhalo, providing a minimum distance from the solar neighborhhod to such a subhalo. So long as the boost factor for the main halo remains sufficiently small, this scenario can therefore overcome the ICS constraints that restricted the standard MH-only model.

### iv.3 Astrophysical prediction and extragalactic constraints

In figure 10 we provide an example of the gamma ray flux predicted by the close subhalo scenario, as compared to the main halo scenario. The gamma ray flux comes predominantly from final state radiation rather than inverse Compton scattering of the annihilation products. We chose the energy bin E = 23 GeV, which is the most constraining for the main halo case. Although both scenarios converge at high latitudes, low latitude measurements have already ruled out the main halo scenario, and provide a way to constrain the model. With more exposure and precise removal of point sources, the Fermi LAT may provide a diffuse background low enough to rule out these predictions. As a further test, census experiments such as the upcoming Gaia satellite may provide a precise enough map of the local gravitational potential to confirm or rule out the presence of such a DM overdensity Brown:2008sh (). Direct measurement of such an overdensity would however be difficult: a subhalo such as SH5, located at a distance that would not saturate gamma ray bounds, would contribute less than of the local DM density.

From previous works, we infer that extragalactic bounds on this scenario are not as strong as the ones we have computed above. Bounds from dwarf spheroidal galaxies could plausibly be important since the velocity dispersions are of the same order as what is required for our subhalo enhancement, i.e. km s Walker:2007ju (). However, the most stringent Fermi LAT bounds Abdo:2010ex () from such galaxies put the upper limit on DM annihilation into a 2 final state at around if only final-state radiation is considered, and around 300 if ICS bounds are included as well. Abdo:2010dk () computed the cosmological dark matter annihilation bounds for the same 2 final state scenario, and find that larger than 300 is excluded at the 90% confidence level. This is using the results of the Millennium II structure formation simulation, and is indeed model-dependent. Extrapolation to the 4 scenario is independent of astrophysics. We can therefore take the results of Meade (); Papucci:2009gd () who have construced bounds on both channels. They show that FSR bounds are consistently an order of magnitude weaker in the 4 case, given the softer photon spectrum in this scenario. We can therefore take these extragalactic results to be far less constraining than the stringent bounds from the center of our own galaxy.

Finally, we verify that this model does not saturate bounds on dipole anisotropy of the cosmic ray spectrum. The dipole anisotropy can be defined as

(21) |

where is the standard dipole power of the measured electron and positron flux in the sky. The Fermi LAT collaboration Ackermann:2010ip () have presented upper bounds on this quantity. These range from at GeV up to at GeV. Given a diffusive model, this can be computed Ackermann:2010ip ():

(22) |

where is the diffusion coefficient (3) and is the density of cosmic ray electrons and positrons, including astrophysical backgrounds. Taking the background to be isotropic, we computed the dipole anisotropy in the case of a single close subhalo producing enough electrons to explain the Fermi excess. In every case falls well below bounds. Results for SH5 are presented in Figure 11. The anisotropy rises monotonically with energy, from 60 GeV (red line) to 500 GeV (black line).

## V Particle physics realizations

In the previous sections we have identified scenarios where subhalos could provide the observed excess PAMELA and Fermi leptons, from a purely phenomenological perspective. In particular, certain values for the annihilation cross section boost factors are needed for the subhalos, and upper bounds for that of the main halo (depending upon assumptions about its density profile) were derived. It is interesting to ask whether simple particle physics models with boost factors from Sommerfeld enhancement can be consistent with these requirements.

The simplest possibility for model building is dark matter that annihilates into light scalar or vector bosons, which subsequently decay into leptons. This class of models automatically gives a boost factor to the annihilation cross section, through multiple exchange of the boson, resulting in Sommerfeld enhancement. However it is not obvious that one can find models with the desired boost factors for the subhalos and main halo. One constraint that limits our freedom is to not exceed the measured density of dark matter. It will turn out that our mechanism works most naturally if the DM responsible for signals in the galaxy is a subdominant component comprising some fraction of the total DM population Cirelli:2010nh (), with .

We focus on the case of a GeV-scale U(1) vector boson that kinetically mixes with the photon. Such models have the advantage of naturally explaining the coupling to light leptons, without producing excess antiprotons that would contradict PAMELA observations. Let us denote the vector’s mass by and the coupling by , with . If is the DM mass, then the Sommerfeld boost factor is controlled by two dimensionless parameters: and , where is the DM velocity in the center of mass frame. A reasonably accurate approximation to the exact Sommerfeld enhancement is given by the expression Cassel:2009wt (); Slatyer:2009vg ()

(23) |

where and . (The cosine becomes if the square root becomes imaginary.)

To take into account leptophilic DM that is only a subdominant component of the total DM, suppose that is the value of that would give the correct thermal abundance, which scales like the inverse annihilation cross section ; then we can parametrize . The rate of annihilations goes like if stands for the leptophilic component of the DM. We accordingly define an effective boost factor

(24) |

where is the intrinsic Sommerfeld enhancement factor. Thus any constraint on in a theory with becomes a constraint on in the more general situation.

### v.1 Averaging of boost factor

Of course, the DM velocity has no definite value; instead we need to average over the possible values within the subhalos and the main halo, weighted by the appropriate distribution function. We take it to be Maxwell-Boltzmann with a cutoff at some escape velocity,

(25) |

This isotropic form is only an approximation since the true
distribution has some small anisotropy between the radial and angular
components; we will for simplicity ignore this complication.
The velocity dispersion depends upon
the radial distance
from the center of the halo or subhalo. The dependence has been
measured for the subhalos in the Via Lactea II simulation; see figure
12. The shape is universal, but is scaled along the
respective axes by parameters and that depend
upon the subhalo. The latter is related to the scale radius
by ; the former is given by (20)
and also listed in table
2 for the subhalos of interest. For numerical
purposes we fit the sides of the curve passing through the
points of fig. 12 by lines (omitting the rightmost point),
and the middle by an inverted parabola.^{5}^{5}5The velocity
dispersion curve is fit by
(26)
where and .
We use the same form of
for the main halo, with kpc and km/s. Other authors have advocated higher values of the
velocity dispersion, km/s at
Savage:2009mk (), which would correspond to km/s in the present parametrization. We will also consider the
higher value to take account of this uncertainty.

The escape velocity can be computed explicitly for the subhalos from the standard result , where is the mass within radius . The result for an Einasto profile is

(27) | |||||

where is the upper incomplete gamma function. For the main halo, this procedure would not be correct because of the significant contribution of baryons, not included here. We adopt the result for of ref. Cirelli:2010nh () for the main halo (see appendix C of that reference).

With these ingredients, we can compute an average Sommerfeld enhancement factor for each subhalo:

(28) |

The factor of in the argument of occurs because the appearing in eq. (23) through is half of the relative velocity. is the appropriate weighting factor because the rate of annihilations is proportional to . For the subhalos, the range of integration for is from 0 to , but for the main halo we take lower and upper limits that correspond to the angular region of the sky that is used to set the gamma ray constraints: kpc and kpc. The reason is that the bound for the main halo comes from the gamma ray constraint rather than from lepton production. We are thus interested in the boost factor relevant to the region of galactic latitude. The distances of closest approach to the galactic center, hence largest rate of ray production associated with these lines of sight, are given by .

### v.2 Relic Density Constraint

The enhancement factor (23) depends rather strongly on the gauge coupling ; therefore it is interesting to know what constraint the relic density places upon . The effect of a Sommerfeld-enhanced DM model on the relic densitie has been discussed by Zavala:2009mi (). Notice that DM transforming under a U(1) gauge symmetry as we have assumed must be Dirac and therefore could have a relic density through its asymmetry, similar to baryons. However, unless the DM was never in thermal equilibrium, then should not be less than the usual value leading to the correct relic density, since otherwise the thermal component will be too large.

There are two kinds of final states for annihilation of DM in this class of models: into a pair of gauge bosons , by virtual DM exchange in the and channels, or into dark Higgs bosons , by exchange of a gauge boson in the channel. Assuming the DM () is much heavier than the final states, the respective squared amplitudes, averaged over initial and summed over final spins, are

(29) |

where is the U(1) charge of relative to (replace for multiple Higgs bosons), is the scattering angle, and we have included the leading dependence on the initial velocity in the center of mass frame. The factor averages to in the integral over . In computing the associated cross section, it must be remembered that the final state consists of identical particles, while the Higgs channel does not. The total amplitude can therefore be written in the form , with

(30) |

if we use the phase space for identical particles.

To find the cross section relevant during freeze-out in the early universe, we thermally average the -dependent following ref. Cline:2010kv (). We include approximately the effect of Sommerfeld enhancement as described there, to obtain

(31) | |||||

The terms that are subleading in , but enhanced by , are due to the Sommerfeld correction. We approximate the freezeout temperature as , the usual result of solving the Boltzmann equation for DM in the TeV mass range, and equate to the value cm/s usually assumed to give the correct relic density. This gives an implicit equation for of the form , which however quickly converges by numerically iterating. Fig. 13 displays the resulting dependence of on for several values of .

The bound that the density of the leptophilic DM component not exceed the total DM density is . We parametrize the coupling by with in what follows.

### v.3 Interpolation between and mixed final states

In our numerical computations with GALPROP, we considered two cases for the final state annihilation channels: either , applicable for gauge bosons with mass , or to a mixture of electrons, muons and charged pions, appropriate for decays of gauge bosons with mass greater than . The relative abundance of , and in the mixed final state can be computed from the branching fractions of the decays, discussed in connection with eq. (8).

For intermediate values of the gauge boson mass, , we can use the branching ratios to interpolate between our maximum-allowed MH or best-fit SH boost factors for the case and those of the fiducial case. The maximum allowed boost factors of the main halo complying with the ICS constraints are taken from table II. To estimate the best fit boost factors for the subhalos in the fiducial final state, we rescale the results shown in table II by the ratio of best-fit boost factors for the main halo, in the MH-only scenario. These ratios are for the Einasto profile and for the isothermal, quite close to unity, and so the best-fit values of the SH boost factors hardly depend upon this scaling. More significant is the change in the best-fit mass, from to 1.2 TeV, which enters into the computation of the Sommerfeld enhancement and the value of the gauge coupling (). We use the branching ratios to interpolate as well. For the MH upper bounds in the small- and large- regions, we use the values from Table II, and interpolate similarly for intermediate .

### v.4 Theoretical fits

For a given value of the gauge coupling , we can determine the predicted boost factors as well as the desired values for each subhalo, as a function of the gauge boson mass , and similarly for the main halo, except here we have an upper bound on rather than a best-fit value. This bound in fact presents the biggest challenge to finding a working particle physics model. For close to the thermal relic density value , the predicted boost factor of the main halo far exceeds the bound , even if we try to decrease by reducing via a large hidden Higgs content or by increasing the dispersion of the main halo. Fig. 14 illustrates the discrepancy for and km/s. Lower values of or only make this tension worse.

As we mentioned above, even though it is not theoretically possible to make the gauge coupling weak enough to solve this problem, ironically one can rescue the scenario by increasing beyond the thermal value, since this suppresses the relic density of the DM component we are interested in, and thus reduces the scattering rate. Allowing decreases both the density of the leptophilic component and the effective boost factor by . With , depending upon the shape of the main halo DM density profile, we can satisfy the constraint on the MH and still have a large enough boost in certain hypothetical nearby subhalos for them to supply the observed lepton excess. The minimum value of that is needed is larger for a cuspy main halo.

We give two working examples in figure 15, one with and km/s (the larger value advocated in ref. Savage:2009mk ()) and assuming an Einasto profile for the main halo, and the other having and km/s (the more standard assumption for the velocity dispersion), with an isothermal halo. In these figures the averaged boost factor of the relevant subhalos are plotted as solid lines, while the required values of are the dashed curves. Wherever these intersect represents a possible value of the gauge boson mass to consistently explain the observed lepton excess. At the same time, the main halo boost factor (lowest solid curve in the small- region) must lie below the black dashed lines to satifsy gamma ray constraints. The rationale for taking the larger value of for the Einasto profile is that larger velocities help to suppress the boost factors and thus make it easier to satisfy the ICS constraint, so that we are not forced to choose an even larger value of . The isothermal profile is less constrained.

In the first panel of fig. 15 with the Einasto profile, only subhalos SH4 and SH5 have large enough boost factors to ever reach the required values. There are many points of intersection, but mainly those for SH5 and in the mass range MeV are consistent with the gamma ray bounds on the main halo. For the isothermal halo, these constraints are less stringent, and it is possible to find points of intersection using for all five of the sample subhalos, although they are much more rare for SH1SH3 than for SH4 and SH5. In this example, the intersection points that respect the ICS bound are restricted to GeV. For larger values of , all the boost factors will be further suppressed, and GeV will become allowed for SH4 and SH5.

One advantage of requiring large is that the corresponding dilution of the DM density by insures that the model satisfies stringent CMB constraints from annihilations in the early universe changing the optical depth Huetsi:2009ex (); Cirelli:2009bb (); Slatyer:2009yq (), as pointed out in Cirelli:2010nh (). The CMB constraint is shown in fig. 16, along with the PAMELA/Fermi allowed regions from ref. Papucci:2009gd () for and final states. The case is allowed by the CMB constraint, but is ruled out. Because our model has at most a fraction of of muons in the final state, it is probably already safe, but the additional weakening of the bound by the factor ensures that this will be the case. Similarly, our scenario overcomes the no-go result of ref. Feng:2010zp (), which pointed out that Sommerfeld enhanced annihilation in the early universe leads to constraints on the MH boost factor which are lower than those needed to explain the lepton anomalies. Our MH boost factor can satisfy these constraints since the MH is no longer considered to be the source of the excess leptons.

The Sommerfeld enhancement is nearly saturated for the low velocities of the subhalos at these large values of , so their curves are nearly overlapping except at the smallest gauge boson masses. The main halo boost factor is not saturated on the other hand, and lies below the FSR bound for most values of . We have chosen the gauge couplings, parametrized by , to nearly saturate the FSR bound. By taking larger (larger ), the bounds could be satisfied by a larger margin. But this would also reduce the values of the subhalos by a similar amount, making it more difficult to get a large enough lepton signal from SH1SH3. SH4 and SH5 would remain robust possible explanations.

## Vi Discussion and Conclusions

We have shown that gamma ray constraints on leptophilic annihilating dark matter are significantly stronger than in previous studies, when we take into account the contributions to inverse Compton scattering from primary and secondary electrons and positrons, before including excess leptons from the DM annihilation. We attribute part of this difference to the method of solving the diffusion equation (1) — fully numerical rather than semi-analytic — meaning that the space-dependence of the diffusion coefficient is taken into account. The difference between the predicted and observed spectra of gamma rays is greatly reduced, leaving little room for new contributions. Because of this, even cored halos, which were allowed by other analyses, become excluded. However, we find that these constraints can be weakened and possibly overcome if annihilations in a nearby subhalo are the dominant source of anomalous leptons, rather than annihilations in the main galactic DM halo. In this way, the PAMELA/Fermi cosmic ray excesses can be explained, without violating bounds from the recent Fermi LAT diffuse gamma ray survey.

It must be admitted that the subhalo loophole we present is rather special. First, only atypically dense subhalos, relative to the sample provided by Via Lactea II, give a large enough boost factor (see fig. 8). Second, the subhalo would need to accidentally line up nearly with the galactic center in order for the ICS gamma rays associated with these leptons to be sufficiently hidden by the noisy background of the GC. Of course, had we neglected ICS contributions of background electrons, similar to previous studies, less fine tuning of the subhalo properties would be necessary. Also we do not require the subhalo to be particularly close; fig. 9 shows that the lepton flux only starts to fall at distances of kpc. Our finding could be regarded as a proof of concept. It is possible that the effects of unresolved substructure within the subhalo Bovy (), which can increase the boost factor, would also make the scenario work more easily. On the positive side, there is the opportunity of testing whether there is such a nearby subhalo, since we predict the spectrum of ICS gamma rays it contributes (see fig. 10). A better understanding of backgrounds, for example from point sources, could make it possible to rule out the proposal.

On the particle physics side, we have shown in detail that the subhalo scenario can be made consistent with one of the simplest models of leptophilic dark matter, where the DM is in a hidden sector that communicates with the standard model only through kinetic mixing with hypercharge of a new gauge boson in the GeV mass range. The relative couplings to leptons and charged pions are completely specified and the model has only two free parameters, the gauge coupling and gauge boson mass (the DM mass is fixed by fitting to the spectrum of anomalous ). The gauge coupling is constrained by the relic density of the DM. The Sommerfeld enhancement factor is completely fixed by and the kinematical halo properties. We find (similarly to ref. Cirelli:2010nh ()) that the predicted boost factor for the main halo is always too large to satisfy ICS constraints unless the leptophilic component of the DM is small, comprising a fraction of order of the total DM. The small fraction can be achieved by assuming is larger than the value required for the usual thermal abundance by the factor . This raises the interesting possibility that the DM that may be responsible for the cosmic ray anomalies is distinct from the dominant DM species that might be discovered by direct detection.

## Acknowledgments

We would like to thank Ilias Cholis for kindly giving us access to his modifications to the GALPROP code, and Troy Porter for valuable help with the Fermi data analysis. We thank Andrew Frey, Thomas Konstandin, Guy Moore, Natalia Toro and Mike Trott for helpful discussions. Our research is supported by NSERC (Canada) and FQRNT (Québec).

## References

- (1) N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, Phys. Rev. D79, 015014 (2009), 0810.0713.
- (2) M. Pospelov and A. Ritz, Phys. Lett. B671, 391 (2009), 0810.1502.
- (3) I. Cholis, D. P. Finkbeiner, L. Goodenough, and N. Weiner, (2008), 0810.5344.
- (4) I. Cholis, G. Dobler, D. P. Finkbeiner, L. Goodenough, and N. Weiner, (2008), 0811.3641.
- (5) P. Meade, M. Papucci, A. Strumia, and T. Volansky, (2009), 0905.0480.
- (6) F. Chen, J. M. Cline, and A. R. Frey, (2009), 0907.4746.
- (7) G. Kane, R. Lu, and S. Watson, (2009), 0906.4765.
- (8) M. Kuhlen, J. Diemand, and P. Madau, (2008), 0805.4416.
- (9) M. Cirelli, M. Kadastik, M. Raidal, and A. Strumia, Nucl. Phys. B813, 1 (2009), 0809.2409.
- (10) J. L. Feng, M. Kaplinghat, and H.-B. Yu, (2010), 1005.4678.
- (11) J. M. Cline, A. C. Vincent, and W. Xue, Phys. Rev. D81, 083512 (2010), 1001.5399.
- (12) M. Cirelli, P. Panci, and P. D. Serpico, (2009), 0912.0663.
- (13) A. Strumia, Prog. Theor. Phys. Suppl. 180, 128 (2010).
- (14) M. Cirelli and J. M. Cline, (2010), 1005.1779.
- (15) J. Bovy, (2009), 0903.0413.
- (16) M. Kuhlen, (2009), 0906.1822.
- (17) M. D. Kistler and J. M. Siegal-Gaskins, (2009), 0909.0519.
- (18) S. Ando, Phys. Rev. D80, 023520 (2009), 0903.4685.
- (19) M. Kuhlen, P. Madau, and J. Silk, (2009), 0907.0005.
- (20) G. Hutsi, A. Hektor, and M. Raidal, JCAP 1007, 008 (2010), 1004.2036.
- (21) M. Kuhlen and D. Malyshev, Phys. Rev. D79, 123517 (2009), 0904.3378.
- (22) P. Brun, T. Delahaye, J. Diemand, S. Profumo, and P. Salati, Phys. Rev. D80, 035023 (2009), 0904.0812.
- (23) A. W. Strong and I. V. Moskalenko, Astrophys. J. 509, 212 (1998), astro-ph/9807150.
- (24) M. Simet and D. Hooper, JCAP 0908, 003 (2009), 0904.2398.
- (25) galprop.stanford.edu.
- (26) T. A. Porter and A. W. Strong, (2005), astro-ph/0507119.
- (27) T. Delahaye et al., Astron. Astrophys. 501, 821 (2009), 0809.5268.
- (28) J. Diemand et al., Nature 454, 735 (2008), 0805.1244.
- (29) V. Springel et al., Mon. Not. Roy. Astron. Soc. 391, 1685 (2008), 0809.0898.
- (30) R. Catena and P. Ullio, (2009), 0907.0018.
- (31) P. Salucci, (2010), 1008.4344.
- (32) P. Salucci and A. Burkert, Astrophys. J. 537, L9 (2000), astro-ph/0004397.
- (33) M. Papucci and A. Strumia, (2009), 0912.0742.
- (34) I. Cholis, L. Goodenough, and N. Weiner, Phys. Rev. D79, 123505 (2009), 0802.2922.
- (35) B. Anderson, M. Kuhlen, R. Johnson, P. Madau, and J. Diemand, Astrophys. J. 718, 899 (2010), 1006.1628.
- (36) A. Birkedal, K. T. Matchev, M. Perelstein, and A. Spray, (2005), hep-ph/0507194.
- (37) G. R. Blumenthal and R. J. Gould, Rev. Mod. Phys. 42, 237 (1970).
- (38) http://fermi.gsfc.nasa.gov/cgi bin/ssc/LAT/WeeklyFiles.cgi.
- (39) T. A. Porter and f. t. F. L. Collaboration, (2009), 0907.0294.
- (40) T. A. Porter, I. V. Moskalenko, A. W. Strong, E. Orlando, and L. Bouchet, Astrophys. J. 682, 400 (2008), 0804.1774.
- (41) I. Cholis et al., (2009), 0907.3953.
- (42) E. A. Baltz and J. Edsjo, Phys. Rev. D59, 023511 (1998), astro-ph/9808243.
- (43) A. G. A. Brown, AIP Conf. Proc. 1082, 209 (2008), 0810.5437.
- (44) M. G. Walker et al., (2007), 0708.0010.
- (45) A. A. Abdo et al., Astrophys. J. 712, 147 (2010), 1001.4531.
- (46) Fermi-LAT, A. A. Abdo et al., JCAP 1004, 014 (2010), 1002.4415.
- (47) Fermi LAT, M. Ackermann et al., (2010), 1008.5119.
- (48) S. Cassel, J. Phys. G37, 105009 (2010), 0903.5307.
- (49) T. R. Slatyer, JCAP 1002, 028 (2010), 0910.5713.
- (50) C. Savage, K. Freese, P. Gondolo, and D. Spolyar, JCAP 0909, 036 (2009), 0901.2713.
- (51) J. Zavala, M. Vogelsberger, and S. D. M. White, Phys. Rev. D81, 083502 (2010), 0910.5221.
- (52) J. M. Cline, A. R. Frey, and F. Chen, (2010), 1008.1784.
- (53) M. Cirelli, F. Iocco, and P. Panci, JCAP 0910, 009 (2009), 0907.0719.
- (54) T. R. Slatyer, N. Padmanabhan, and D. P. Finkbeiner, Phys. Rev. D80, 043526 (2009), 0906.1197.
- (55) G. Huetsi, A. Hektor, and M. Raidal, (2009), 0906.4550.