Cosmic Ray Anomalies from the MSSM?

The recent positron excess in cosmic rays (CR) observed by the PAMELA satellite may be a signal for dark matter (DM) annihilation. When these measurements are combined with those from FERMI on the total ($e^++e^-$) flux and from PAMELA itself on the $\bar p/p$ ratio, these and other results are difficult to reconcile with traditional models of DM, including the conventional mSUGRA version of Supersymmetry even if boosts as large as $10^{3-4}$ are allowed. In this paper, we combine the results of a previously obtained scan over a more general 19-parameter subspace of the MSSM with a corresponding scan over astrophysical parameters that describe the propagation of CR. We then ascertain whether or not a good fit to this CR data can be obtained with relatively small boost factors while simultaneously satisfying the additional constraints arising from gamma ray data. We find that a specific subclass of MSSM models where the LSP is mostly pure bino and annihilates almost exclusively into $\tau$ pairs comes very close to satisfying these requirements. The lightest $\tilde \tau$ in this set of models is found to be relatively close in mass to the LSP and is in some cases the nLSP. These models lead to a significant improvement in the overall fit to the data by an amount $\Delta \chi^2 \sim 1/$dof in comparison to the best fit without Supersymmetry while employing boosts $\sim 100$. The implications of these models for future experiments are discussed.


Introduction
Evidence for physics beyond the Standard Model (SM) has been uncovered with the cosmological observations that ∼ 25% of our universe is comprised of Dark Matter (DM). Numerous astrophysical observations point to the existence of cold DM that cannot be accommodated in the SM [1]. This cold DM could take many forms, for example a weakly interacting massive particle (WIMP), axion, dark photon, or something more exotic are all compatible with the data [2]. The WIMP-like DM scenario has received the most consideration as WMAP measurements [3] of the relic density point to the TeV-scale for weak-strength DM annihilation and WIMPs occur naturally in many theories which go beyond the SM. As a result, much attention is focused on various means of detecting WIMPs, through direct detection experiments deep underground, direct production at colliders, or indirect detection in cosmic ray (CR) spectra.
In fact, an anomaly has recently been reported in CR data. The PAMELA collaboration observes [4] a steep rise in the positron fraction (the ratio of positrons over the sum of electrons plus positrons) at energies between 10 and 100 GeV. The Fermi Gamma Ray Space Telescope (FERMI) reports [5] a slightly stiffer electron spectrum than expected, but with no dramatic features. These are countered by observations of the anti-proton fraction (ratio ofp to p) by PAMELA [6] and various gamma ray measurements by FERMI [7,8,9,10] which reveal no unexpected characteristics.
Much attention has been given to explaining these results in theories with new physics, with Supersymmetry (SUSY), arguably the most widely-studied theory of physics beyond the SM [11], receiving the most consideration. In these models the DM particle is identified with the Lightest Supersymmetric particle (LSP) which is typically assumed to be the lightest neutralino. However, it was demonstrated early on [12] that the most-studied Supersymmetric models, such as minimal Supergravity (mSUGRA), coupled with a thermal cosmology, have difficulty accommodating this data. In this case, separate issues arise concerning the reproduction of the sharp rise in the positron fraction, while simultaneously not affecting the anti-proton or gamma ray spectra. In addition, the typical Supersymmetric thermal DM annihilation cross section is found to be much too small, by a factor of order 10 3−4 , to explain the PAMELA positron spectrum excess. This requires the introduction of a large boost factor, B, which multiples the annihilation cross section, hence increasing the rate. It has been hypothesized that such a boost factor could arise from clumpiness [13] in the DM distribution, from resonant [14] or Sommerfeld enhancement [15] of the annihilation cross section, from non-minimal SUSY contributions [16] [17], or from a non-thermal cosmology [18]. Meanwhile, once such a boost factor is incorporated it can create problems with the SUSY contributions to both thep/p and gamma ray fluxes, rendering them too large [12].
In this work, we study a broader and more comprehensive region of the parameter space than is usually studied in the Minimal Supersymmetric Standard Model (MSSM) with the expectation of obtaining a variety of contributions to cosmic rays that differ from those found in mSUGRA. We employ the large set of models (∼ 68.5k) generated in the 19dimensional phenomenological MSSM (pMSSM) parameter space by Berger et al. [19]. These models contain no theoretical prejudice concerning the origin of Supersymmetry breaking or interactions at the GUT scale. As noted in this previous work, the characteristics of these pMSSM models can be quite different from those of mSUGRA. In this paper we examine the predictions of these models for the various CR spectra in the hope of obtaining a good fit to the data without the need to incorporate such very large boost factors.
To this end, we couple this scan of the pMSSM parameter space with an exploration of the uncertainties inherent in the calculation of the astrophysical contributions to cosmic rays. There are numerous steps in the computation of the propagation of cosmic rays across the galaxy, each fraught with their own approximations and uncertainties. We account for these by performing an extensive scan over the astrophysical parameters entering the propagation calculation and find a set of ∼ 6k phenomenologically viable CR 'models' that are consistent with the CR data set (except for the high energy bins in the PAMELA positron flux) and are subject to a number of additional constraints. We then compute both the standard astrophysical background spectra and the SUSY DM annihilation contributions using the same propagation parameters, for each combination of the pMSSM and astrophysics models. Performing a global χ 2 fit to the (e + + e − ), e + /(e + + e − ) and (p/p) data, we find a region of pMSSM parameter space that yields a good fit to the positron and anti-proton fractions as well at the (e + + e − ) data (e.g., χ 2 ∼ 1.55/dof) with boost factors in the range ∼ 100 − 200, far lower than those typically required in conventional mSUGRA. The presence of the SUSY DM component is found to significantly increase the quality of the fit to this data, by ∼ 1 unit of χ 2 /dof, in comparison to the null hypothesis of no new physics with only conventional astrophysical sources. We next perform a consistency check, employing data from both the mid-latitude diffuse gamma ray spectrum as well as from the total γ flux from faint dwarf galaxies for a small subset of 10 models which are consistently found to provide the best overall fit. In these models we find that the pMSSM spectrum and couplings are such that the LSPs tend to annihilate almost exclusively into tau pairs. This mostly occurs because: (i) the LSPs are found to be dominantly bino-like, preventing annihilations into W, Z-bosons, (ii) the staus are light, (iii) while the sbottoms are reasonably heavy thus suppressing annihilation into bb and avoiding significant contributions to the antiproton flux.
The rest of this paper proceeds as follows: the next section describes the generation and properties of the pMSSM model sample that was originally carried out in [19]. We discuss our CR results in Section 3, first describing the CR propagation uncertainties and our scan of the astrophysics parameter space, next, giving the results that are obtained from astrophysical contributions alone, and, finally, including the SUSY contributions. In Section 4 we compute the SUSY contributions from WIMP annihilation in dwarf galaxies as well as there contribution to the mid-latitude diffuse gamma ray flux, while Section 5 contains a description of the best-fit pMSSM models. Our conclusions are given in Section 6.

Generation and Properties of pMSSM Model Set
We study a large set of Supersymmetric models that broadly sample the parameter space, do not incorporate a specific SUSY-breaking mechanism and are free of theoretical assumptions at the GUT scale. We do so in the expectation that a more extensive sampling of the MSSM parameter space will be better able to identify models which accommodate the cosmic ray data. We use the MSSM database generated in the work of Berger et al. [19]. In this section, we briefly describe the procedure employed in constructing these models and summarize the properties of these models that are most relevant to DM phenomenology. Throughout this work, we refer to a point in the Supersymmetric parameter space as a model.
The lower bounds were chosen for general consistency with collider data and the upper bounds were taken to ensure large sparticle production cross sections at the LHC. A smaller number of points were scanned using logarithmic priors to ensure that the results displayed little dependency on the choice of priors. These points were then subjected to a set of theoretical and experimental constraints summarized below. The number of points scanned and analyzed was limited by CPU availability.
The theoretical requirements imposed on models were that their SUSY spectra contain no tachyons, that there be no color (or charge) breaking minima in the scalar potential, that electroweak symmetry breaking be consistent and that the Higgs potential be bounded from below. Consistency with the electroweak precision measurements ∆ρ, the invisible width of the Z-boson, and the anomalous magnetic moment of the muon was required. Constraints from the flavor sector were applied and agreement with the rates or limits for b → sγ, B s → µµ, B → τ ν and meson mixing was imposed. Limits from collider searches for sparticles and SUSY Higgs bosons were enforced. The charged sparticle and SUSY Higgs bounds from LEPII contain many caveats and dependencies on the details of the individual model spectrum; these were carefully incorporated into the analysis. The Tevatron searches for squarks and gluinos were generalized to render them model independent and multi-jet + MET events were generated for each pMSSM model and compared to the data. In addition, bounds from the stable charged particle searches at LEPII and the Tevatron were applied. Two astrophysical constraints were imposed on the long-lived relicχ 0 1 . Agreement with the WMAP 5-year measurement of the relic density was required such that Ωh 2 | LSP ≤ 0.121 (1σ above the central value). In not employing a lower bound on Ωh 2 | LSP , the possibility was left open that dark matter may have multiple components within or outside of the pMSSM and thermal relic framework. Lastly, restrictions on spin independent and spin dependent cross sections from direct DM detection searches were incorporated. For further details and references on this set of constraints, see [19].
After this set of constraints was imposed, ∼ 68.5k pMSSM models from the flat prior sample were found to survive, satisfying all the restrictions. Here, we discuss the attributes of these models which are most germane to the present analysis. Figure 1 (reproduced here from [19]) presents a histogram of the masses of the four neutralino species for this model set. We see that the LSP mass lies mostly in the range 100 − 250 GeV in these models. Models with a mostly Higgsino or Wino-like LSP generally have a chargino with nearly the same mass as the LSP, and as sufficiently light charginos would have been detected at LEP or the Tevatron, there are fewer models with an LSP of mass mχ0 1 < ∼ 100 GeV. The gauge eigenstate content of the LSPs in these models is described in Table 1. We emphasize that the percentages given in this Table apply strictly to this pMSSM model sample and are not related to the likelihood of models being realized in Nature. We note that mostχ 0 1 's are relatively pure eigenstates, with models where the LSP is Higgsino or mostly Higgsino being the most common case. Within mSUGRA, the LSP is generally close to being a pure Bino, suggesting that this pMSSM model set has substantially different properties that will affect DM annihilation rates.
The identity of the next-to-lightest Supersymmetric particle (nLSP) is quite varied in this pMSSM model sample. The lightest chargino is the nLSP in roughly 78% of the models. This is because a large fraction of models have a Wino or Higgsino eigenstate LSP (as seen in Table 1) and there is generally a small mass splitting between a mostly Wino/Higgsino neutralino and the corresponding chargino. The second lightest neutralino is the nLSP ∼ 6% of the time. These are generally models where the LSP is dominantly Higgsino. Whileχ ± 1 andχ 0 2 are the nLSP in the vast majority of cases, in the remaining 16% of the model sample 10 other sparticles are the nLSP with approximately equal probability. These are theτ 1 ,ẽ(µ) R ,ν τ ,ν e(µ) ,b 1 ,t 1 ,ũ R,L ,d R andg. Lastly, two models have thed L as the nLSP. Scenarios in which these sparticles are the nLSP may lead to uncharacteristic DM annihilation rates in various channels (such as a large rate into τ pairs), as well as interesting signatures at the LHC [21]. Many models are found to have small mass splittings between Figure 1: Distribution of neutralino masses for the pMSSM model sample. From [19]. content. This is defined by the modulus squared of the elements of the neutralino mixing matrix in the SLHA convention. See [11] for details.

LSP
the LSP and nLSP, leading to potentially large DM annihilation rates.
As stated above, it was not required that the LSP account for all of the dark matter. Rather the only constraint was that the LSP relic density not be so large as to be inconsistent with WMAP. Figure 2 (reproduced here from [19]) displays the distribution for the LSP relic density in our pMSSM model set. Note that this distribution is peaked at rather small values of Ωh 2 | LSP . In particular, the mean value for this quantity is Ωh 2 | LSP ∼ 0.012. There are 1240 models that produce a relic density in the range 0.10 ≤ Ωh 2 | LSP ≤ 0.121, thus saturating the WMAP measurement; we will refer to these hereafter as the high-Ω model set. For many of the observables that we compute here, we will need to scale by the ratio This gives the fraction of the local relic density as determined by WMAP that arises from Supersymmetric DM. In our numerical computations, we take the central value Ωh 2 | W M AP = 0.1143. The importance of including this scale factor is demonstrated in Fig. 3, which shows both the scaled and unscaled total annihilation cross sections for our pMSSM model set. The points highlighted in orange correspond to the high-Ω models. We see from the figure that this scale factor is necessary for consistency with WMAP, which bounds the total annihilation cross section. The quantity σv R 2 is more indicative of the annihilation signals that result from a given model. From [19]. Of particular interest to us below, are the relative DM annihilation rates into τ pairs versus bb final states in our pMSSM model set. This is presented in Figure 4 which shows the ratio of these two final state annihilation rates versus the scaled rate for σv τ + τ − for models with mχ0 1 > 100 GeV. The points highlighted in orange correspond to the high-Ω models. We note that for the full model sample, the values of the ratio σv τ + τ − / σv bb range over more than ten orders of magnitude. We see that there are many models with a high annihilation rate into τ pairs at relatively large rates with high purity (i.e., the rate into bb is significantly lower). These models are perhaps the best candidates for "leptophilic" DM contributions to CR spectra necessary to simultaneously explain the (e + + e − ) andp/p fluxes observed by PAMELA. In the analysis that follows, we will repeatedly find reasons to focus on this subset of models.

Cosmic-Rays
We now turn to a discussion of cosmic-ray (CR) signals from the annihilation of WIMPs in our pMSSM model set. We pay particular attention to the anomalous measurements of CR e ± species (here termed, "CRE"), investigating the properties of models that do best to alleviate the tension between measurements of the (e + + e − ) , e + /(e + + e − ) andp/p flux spectra.
It is generally expected that the usual thermal cosmological evolution of the MSSM particle content should give a DM annihilation cross-section of σv R 2 < ∼ 3 · 10 −26 [22], a number well shy of that required for an easy fit to the PAMELA measurement of the e + /(e + + e − ) spectrum. This number, though necessarily somewhat approximate, has served to crystallize the lore that the CRE anomaly must be described by non-standard DM models (if by DM at all). In this work we test this belief by performing an extensive exploration of both the parameter space describing phenomenologically-viable supersymmetric models and the space of parameters used to model the propagation of cosmic-rays in our galaxy. We will look for regions in both spaces of parameters which provide a good fit to all relevant data.
We will follow the now typical practice of adding a so-called boost factor, that is: σv R 2 → B · σv R 2 . While boost factors have been motivated by e.g., non-thermal cosmological history [18] or novel halo distribution effects [13], here we will remain agnostic as to their origin and just quantify the boost that is necessary to best fit the data with respect to the usual scenario: a thermal cosmology and the canonical value for the local dark matter density (0.3 GeV cm −3 ). We use the same boost factor in adding the SUSY contributions to all measurements, for example in the e + /(e + + e − ) and inp/p fluxes 2 . It is not difficult to find MSSM models for which boosts of O(10 3 − 10 4 ) allow good fits to the data. However, it has been argued that effects in the halo distribution alone can only be expected to account for a boost factor of ∼ 2 − 3, due to the uncertainty in the estimation of the local halo abundance [23]. A possible additional factor < ∼ 10 could arise from substructure within the Milky-Way halo [24], and, for signals from dwarf galaxies, there could possibly be a factor < ∼ 100 [25]. In this work we will focus on models that best fit the data with B < 200.
Of course the size of B, and the effect of SUSY annihilations in general, can only be sensibly determined if one has a good estimate of the relevant astrophysical backgrounds. In this respect, the treatment of CRE spectra is particularly difficult due to the large uncertainty inherent in modeling their sources and propagation [26] [27]. In this section we discuss how we take account of these uncertainties by performing a large scan over some of the astrophysical parameters used to model cosmic-ray propagation, thereby creating a variety of astrophysical background spectra that we regard as internally consistent and plausible with respect to relevant astrophysical measurements.
As mentioned in Section 2 the LSPs in our model set typically have mχ0 1 < ∼ 500 GeV, well shy of what would be required to provide a DM explanation to the putative anomaly seen by FERMI-LAT in the (e + + e − ) spectrum near 1 TeV. As has already been emphasized [28], the FERMI-LAT measurement of the (e + + e − ) flux can be well described in a conventional diffusion model if the injection spectra of CREs coming from the average background of CR sources is stiffer than previously thought, though this interpretation only exacerbates the apparent anomaly in the PAMELA observation of the e + /(e + +e − ) spectrum. Here we search for a region of astrophysical and pMSSM parameter space such that the excess measured by PAMELA in the e + /(e + + e − ) spectrum is alleviated by a sizeable contribution of positrons from pMSSM DM annihilation while at the same time the FERMI-LAT measurement of the (e + + e − ) flux is accommodated by the astrophysical background model (without including any particular astrophysical sources, such as pulsars).
At the end of this section we combine the two scans, over both astrophysical and pMSSM parameter spaces, and perform a simultaneous fit of the PAMELA e + /(e + + e − ) and (p/p) flux measurements as well as the FERMI-LAT (e + + e − ) flux measurement. We then discuss the properties of both the astrophysical models and the SUSY models that lead us to the best global fits with the smallest boost factors.

Consistent Modeling of the High-Energy Sky
The prevailing theory describing CR propagation throughout the galaxy is based on the diffusion model, for which a detailed description can be found in [29] [30]. The diffusive paradigm follows naturally from the expectation that charged CR species follow tangled paths as they scatter resonantly off of the turbulent features in the galactic magnetic field whose size match their gyroradii; this explains the observed isotropy of CRs and the retention of CR species in the galaxy. In its simplest guise one can solve for the local spectra of a single CR species using a model with relatively few parameters. The density of a CR species at galactic radius r and with momentum p, i.e. ψ(r, p, t), can be calculated using the CR transport equation [29]: One solves this equation over a diffusion zone that is usually represented by a thick cylinder whose origin and axis coincide with that of our galactic disk and whose half-height is typically in the range L ∼ 1 − 10 kpc. The diffusive term has an energy-dependent diffusion coefficient, D xx = D xx ( r, p), that describes scattering off of the turbulent features of the galactic magnetic field, δB( r, p) (the full field being thought of as the sum of smooth and turbulent components, B( r, p) = B sm ( r, p)+δB( r, p)). One typically simplifies to a spatially uniform diffusion coefficient with free-escape boundary conditions (i.e. ψ(r, p, t) = 0) at the edges of the diffusion zone. The other terms describe effects such as energy losses by radiation (ṗ), diffusion in momentum space (so-called "re-acceleration," D pp ), interactions of single CR particles with the convective "winds" created by the bulk population of galactic CRs (parameterized by the convection velocity, V c ), and the possible fragmentation and radioactive decay of unstable nuclei (τ f and τ r ), as well as the source term, q( r, p, t). We will describe some of these effects in more detail below.
One can often use the transport equation and simplify the problem enough to obtain explicit analytic solutions for individual CR spectra. However, a network of relations connect various aspects of the astrophysical environment and the spectra of CRs/γ's, making a more holistic analysis desirable. For example, primary protons 3 , the most abundant CR species, source all of the other species of cosmic-rays (including e ± ) by spallative interactions with dust and gas. CRE are strongly related to a vast spectrum of photons in the galaxy as they generate copious amounts of radio photons by synchrotron radiation and gamma-rays by Inverse Compton (IC) up-scattering of ambient photons in the microwave-to-optical range. The parameters in Equation (3) represent a simple phenomenological model that encodes the effects of these relations on a given CR species. However, arbitrary tuning of these parameters without consideration of the effect on related quantities may lead to meaningless results. Using CREs as an example, one may wish to change the term describing electron energy losses in order to see the effect on locally measured CRE spectra while forgetting the consequences of this change for previously measured quantities such as diffuse gamma-ray (via IC scattering) or radio photon spectra (via synchrotron radiation).
The GALPROP numerical package is by far the most developed tool for consistent cosmic-ray analysis [31] [32]. The code solves a network of transport equations for Z≥1 nuclei as well as for electrons and positrons while computing energy losses using realistic maps of galactic gas [33] and of the far-infrared, optical and CMB photons that make up the Inter-Stellar Radiation Field (ISRF) [34]. One also has the ability to transfer tabulated Green's functions from GALPROP to DarkSUSY 5.0.4 [35], where the injected spectra of CREs from dark matter annihilations are calculated using all of the details of the SUSY model spectrum and couplings. In this way we can propagate the CRE coming from dark matter annihilations in the same manner as those arising from standard astrophysical sources so that signal and background may be combined in a consistent manner.

Designing Astrophysical Models
Here we investigate astrophysical background uncertainties by carrying out a scan over the parameters associated with cosmic ray propagation. Our scan is not "without prejudice," in the sense that we do not scan all parameters over their full range of plausible values. In principle there are many combinations of parameters that do not allow for a simultaneous fit of CRE measurements, with or without a SUSY contribution. We instead focus on designing astrophysical models that are both plausible from the sense of astrophysical measurements and that have the potential for a good fit with a supersymmetric DM contribution and a low boost factor. In this Section we describe the process of generating such models.
We use GALPROP av50.1p and in all cases we generate astrophysical models by starting with a conventional model 4 and changing only the parameters that are explicitly mentioned below.
The connection between locally-measured CR spectra and the nature of their interstellar propagation is confounded by the fact that charged particles are heavily affected by the currents and magnetic fields of the heliosphere. This so-called solar modulation effect is more important for low energy CRs (see e.g., [37]) and its modeling is currently an active area of CR research. In this work, the inter-stellar spectra of all species are solar modulated according to a spherical force field approximation [38] with potential φ = 500GV, though in all cases we fit only to data in the range E > 10 GeV to minimize the importance of uncertainties in modeling this effect. This method of computing modulation effects does not depend on the sign of charged species, and so factors out in ratios such as e + /(e + + e − ) , but nevertheless has some effect on fitting absolute spectra (e.g., (e + + e − ) ). We note that, while we minimize the effect of solar modulation on our fit to the positron fraction by simply leaving out the data bins below 10 GeV, previous works (for example, [16]) have attempted to model charge-sign-dependent solar modulation using the data presented in [39]. We make no attempt at charge-sign-dependent solar modulation in this analysis but note that recent work [40], using data that is significantly more recent that that of [39], has shown that the various positron fraction data sets, as well as the predictions of diffusion models such as GALPROP, can be brought into accord (at energies below 10 GeV) in this way.
The astrophysical parameters that we scan can be classified into related subsets: those describing nucleonic CR sources, those that model the diffusion of CRs, those describing electron CR sources and those that model electron energy losses. In each of these subsets one is constrained by associated measured quantities: the proton absolute flux andp/p spectrum are directly related to the nucleonic CR source distribution, B/C and 10 Be/ 9 Be flux spectra ratios are very sensitive to changes in the parameters describing diffusion, and the (e + +e − ) , e + /(e + +e − ) and diffuse mid-latitude gamma-ray spectra are all strongly affected by changes in the electron source and energy loss models. We will describe our treatment of each of these subsets of parameters in the subsections that follow. The phenomenological fact that Z≥1 nuclei undergo negligible radiative energy losses leads to a decoupling of the hadronic quantities (e.g.,p/p , B/C and similar ratios) from the astrophysical parameters that affect only the CRE spectra (e.g., (e + + e − ) , e + /(e + + e − ) and similar quantities), so there is a logical sequence for scanning parameters using minimal CPU time. This approach proceeds as follows.

Nucleon Source and Turbulent Diffusion Parameters
As protons are by far the most numerous CR species, their absolute spectrum has been observed by many experiments, producing data with very little statistical error. We thus begin by finding the values of the nucleon source parameters that best fit this data. We use the GALPROP default spatial distribution for nucleon sources and tune the overall normalization and spectral index (i .e., the exponent of the power-law describing the energydependence of the sources) of these sources to best fit measurements of the proton absolute spectrum from BESS/BESS-TeV [37] [41], AMS-01 [42], CAPRICE98 [43] and ATIC [44]. It is clear that the data favor a (locally-measured) spectral index of ≈ 2.75 ± 0.05, though the ATIC points stiffen somewhat above ∼ 1 TeV; including all four data sets, we find a best fit spectral index of γ n = 2.73 and normalization 3.96 × 10 −6 GeV −1 cm −2 sr −1 s −1 @ 100 GeV. Once the normalization and spectral index are fit in this way, these values are held fixed for the rest of our analysis. The diffusion sector is described by five parameters: D 0xx , δ, L, V a and ∂V c /∂z. Here, D 0xx and δ parameterize the energy-dependent diffusion coefficient as 5 is taken to be spatially constant inside the diffusion zone) and L is the half-height of the diffusion zone. V c describes the convection of single CRs caused by the bulk population of CRs and is taken to be directed vertically away from the galactic disk with constant gradient, ∂V c /∂z. V a is the Alfven velocity; it is related to D pp from Eq. (3) via the approximate relation D pp = p 2 V 2 a /(9D xx ), and so parameterizes the effect of diffusion in momentum space.
We can gain some intuition for the effects of these parameters in the context of a simple 1D diffusion model. We can simplify Eq. (3) by specializing to stable nuclei (τ r , τ f → ∞ andṗ ≈ 0, as nuclei undergo negligible radiative energy losses) and dropping the terms involving V a and V c , which have a small effect on spectra above ≈ 20 GeV [26] [27]. We are then left with only the spatial diffusion term. Taking the source q(z, E) = N δ(z)E −γ , i .e., a source localized to the galactic plane having spectral index γ and normalization N , Eq. (3) becomes: With free-escape boundary conditions, ψ(±L, E) = 0, this is easily solved. We may also obtain a solution for secondary CRs (produced by spallation events involving primary CRs) by substituting the primary spectra as a source for secondaries. At our location (z = 0) we have where the source spectral index, γ = γ n , is universal in the GALPROP model for nuclei of all (Z,A). Secondary species are injected with index γ n + δ (as they are formed by proton spallation events) and σ represents the emissivity for primaries moving through the disk. Although this result was obtained in a somewhat approximate fashion, it is very close to the numerical results obtained using GALPROP to propagate stable nuclei. From this we note that the ratio of secondary to primary spectra is expected to scale as (secondary)/(primary) ∼ (L/D 0xx )E −δ , and thus, for example, the ratio of the measurements of boron CRs (which are known to be almost entirely secondarily produced) to that for carbon CRs (which are known to be dominantly primary in origin) may constrain plausible values of δ and the ratio (L/D 0xx ). In fact the B/C flux ratio is quite sensitive to the choice of diffusion parameters and can tightly constrain candidate diffusion models, as it has been measured in many experiments with good statistics. We describe our scan in this space and the results of our fit to the B/C flux data below, but first, it is worthwhile to make one more observation about this simplified 1D result.
To consider the propagation of CRE species, which lose energy rapidly through radiative processes (e.g., synchrotron radiation and IC scattering), the 1D model above should also include the radiative energy loss term. Despite this, we can get a feel for how to choose δ for our purposes here, by examining this diffusion-only result. To the extent to which the flux ratio e + /(e + + e − ) ∼ e + /e − (CRE electrons far outnumbering CRE positrons), the background-only positron ratio scales as e + /(e + + e − ) ∼ E −γn+γe−δ . Here, the primary electrons are injected with source spectral index γ = γ e , while the positrons (being formed predominantly by proton CR spallation) are injected with spectral index γ = γ n + δ. To this background we would like to add a primary positron contribution from SUSY DM annihilation; this contribution will scale as e + susy ∼ E −γsusy−δ . Hence, from both the standpoint of attaining large e + fluxes from DM annihilations, and from the standpoint of generating e + /(e + + e − ) flux backgrounds that will not require large boost factors when combined with the DM signal, we would like to choose δ to be as small as possible.
As δ describes the spectrum of turbulent features in the bulk magnetic field, there has been a large industry devoted to theoretical and experimental considerations of what its value should be [45] [46]. Particularly well-motivated theoretical values include δ = 0.33 (Kolmogorov turbulence) and δ = 0.5 (Kraichnan turbulence). Experimentally, there are measurements of B/C and various phenomenological propagation models that are used to constrain δ given the measured data. These works have found ranges of values for δ that provide consistency with B/C measurements, some favoring lower values δ = 0.3 − 0.5 [45] and others finding that larger values δ = 0.5 − 0.8 are preferred [46]. In this work we use only propagation models that have δ = 0.33.
One should also note that it is not entirely trivial to identify the "δ" that is being measured by the B/C data with the δ parameter that effectively describes the diffusion of CRE species. This is because CRE energy losses are dominantly radiative, with rates proportional to the Thomson cross-section σ T , whereas the corresponding radiative loss rate for CR nuclei is proportional to (Z 2 /A) 2 (m e /m p ) 2 σ T . Such losses are thus completely negligible for nuclei such as boron and carbon, which lose energy predominantly via interactions with matter, and so should be able to diffuse substantially throughout the galactic magnetic field from their point of creation. Conversely, high-energy ( > ∼ 10 GeV) e ± lose energy very efficiently through synchrotron radiation and IC scattering and so must come from a region within the local ∼ kpc to be detected near our solar system [26] [27]. In such a picture the "δ" appropriate to fitting B/C effectively describes an average over the entire diffusion zone while the δ parameter appropriate for modeling CRE flux spectra near earth is that which describes diffusion in the local galactic neighborhood, a scale at which diffusion is likely highly complicated 6 .
We fix δ = 0.33 and scan over the remaining diffusion parameters, D 0xx , L, V a and ∂V c /∂z, to find configurations that best fit the B/C data. Although there are many measurements of B/C spectra at various energies (see e.g., the list presented in [48]), we choose to fit only the data sets from the CREAM [49], HEAO-3 [50] and ATIC experiments [51]. These sets are the most recent, cover the largest range of energies and have the smallest quoted errors 7 .  Table 2: Best fit parameter configurations in the diffusion sector.
As noted above, B/C provides quite a sensitive constraint on the ratio (L/D 0xx ), but cannot separately constrain either L or D 0xx . This degeneracy can be broken by measuring the abundance of unstable radioactive nuclei relative to their stable isotopic counterparts, i .e., flux ratios such as 10 Be/ 9 Be, which are referred to as "radio-clocks." There are a relatively small number of measurements of such ratios, however. Rather than trying to fit 10 Be/ 9 Be directly in this analysis, we note that the data seem to favor L >  [29]). For each value of L, we scan over the remaining parameters, D 0xx , V a and ∂V c /∂z, and compute a χ 2 fit to the three data sets for each configuration (taking only the bins which cover energies above 10 GeV to avoid the effects of solar modulation), first using a coarse grid and then zeroing in with a finer scan. The resulting best fit configurations are listed in Table 2.
It should be noted that ignoring data below 10 GeV, as we do here, likely has an impact on the determination of best fit parameters, especially those describing convection and re-acceleration. However, this should have little impact on our analysis as the effects of convection and re-acceleration are largely unimportant to CRE spectra above 10 GeV, at least for commonly accepted sizes of these effects [26] [27]. Our analysis is similar in this regard to that of [48], especially in the finding that ∂V c /∂z ≈ 0. We note that one should refer to more serious attempts at fitting B/C in the ∼ 1 GeV range for applications that may be particularly sensitive to these effects.

Electron Source and Energy Loss Parameters
In the case of nucleonic CRs, where measurements of protons have been made with large statistics, and for which propagation modeling is relatively simple, the spectral index of the sources can be relatively unambiguously connected to the locally measured spectral index 8 . In the case of CRE, however, energy losses play an important role in propagation and there is much greater uncertainty in relating locally measured fluxes to the properties of the sources. We vary the nature of the electron sources in a manner similar to that for the nucleon sources, scanning only the overall normalization and spectral index. The normalization of the primary electron sources is treated somewhat differently than in the nucleon case in that we allow an a posteriori adjustment of the normalization of the source to best-fit the data, as we will describe below. We study the nature of CRE energy losses by scanning over the overall rates of synchrotron and IC radiative losses by tuning the galactic magnetic field and inter-stellar photon densities, respectively. We vary the nature of IC energy losses, i .e., whether the IC scatterings are dominantly in the Thomson regime or dominantly in the Klein-Nishina regime, by changing the average energy of inter-stellar photons (in practice, the relative amounts of far-infrared versus optical photons). In this section we give a detailed description of this process.
Since primary electrons dominate over other (secondarily produced) CRE species we expect that the fluxes behave as (e + +e − ) ≈ e − and e + /(e + +e − ) ≈ e + /e − , at least at energies > ∼ 10 GeV. Hence, the choice of primary electron injection index affects each quantity in its own way: decreasing γ e will stiffen the combined (e + + e − ) spectrum while simultaneously softening the energy dependence of the positron fraction. As has been noted previously [26], the uncertainty inherent in γ e can drastically affect the visibility of a potential DM annihilation signal in the positron fraction. The FERMI-LAT measurement of the local (e + + e − ) spectrum is consistent with a power-law of spectral index γ local ≈ 3.05. As was shown in [28], this result can be interpreted in the context of a simple diffusion model with a relatively stiff γ e = 2.42 ("Model 1" of [28]). In such a scenario, however, the expected astrophysical contribution to the positron fraction is even softer than what was previously predicted, diving away from the PAMELA e + /(e + + e − ) flux data, and further exacerbating the anomaly.
Here we scan CRE source spectral indices in the range: 2.42 ≤ γ e ≤ 2.60. This range is chosen because we want astrophysical backgrounds that will require relatively low boost factors when including a DM signal. We will rely on the scan of parameters in the radiative loss sector to make up for the softer spectrum arising from this range of γ e values to obtain the measured dependence, (e + + e − ) ≈ E −3.05 .
We also vary the overall normalization of the CRE source, N e . However, we do not do this from the outset (i .e., in the galdef input files used in GALPROP) as we do for the other parameters, but rather during the fitting process and in the calculation of the χ 2 fits. This is possible as the overall normalization factors out of the numerical process in GALPROP. Thus N e can be treated in the same manner as the SUSY boost factor: N e is allowed to float in determining the best-fit χ 2 and we subtract an additional degree of freedom in calculating the reduced χ 2 .
Let us now discuss CRE energy losses in more detail. At the energies with which we are concerned (∼ 10 − 100 GeV) CRE losses dominantly occur through IC up-scattering of ambient photons and through the synchrotron radiation induced in transiting the large scale galactic magnetic field. We can write a general expression for such losses as [54]: where the first term describes losses via synchrotron radiation, witḣ where E = γm e c 2 and u B = B 2 /8π. The second term describes IC energy losses, with u 0 being the distribution of ambient photons, 0 is the ambient photon energy scaled as 0 = hν 0 /m e c 2 and u 0 = u 0 d 0 . This term would have the Thomson-like energy dependence present in the synchrotron term (i .e.,Ė syn ∼ E 2 ), if not for the energy dependence of the kernel f KN (E, 0 ). An analytical solution is known for f KN (E, 0 ) [55], but for our purposes here we merely summarize some key features: In the idealized case of a monochromatic photon ( 0 ) bath and a low-enough magnetic field, this transition of f KN (E, 0 ) would become manifest in the CRE spectra as a "step" feature [56]. For CREs with energy high enough (w.r.t. photons of energy 0 ) to be in the KN regime, the spectral index stiffens relative to the Thomson result. As energies increase, the overall energy loss rate goes from being dominantly due to IC scattering to being dominantly due to synchrotron radiation (whose contribution always grows asĖ syn ∼ E 2 ), and the CRE spectral index returns to a Thomson-like slope (with normalization larger by a factor u 0 /u B ). The transition between these two regimes occurs at 0 ≈ (1/4γ), corresponding to photons of energy ≈ 1 eV (≈ 10 eV) for electrons of energy ≈ 100 GeV (≈ 10 GeV), respectively.
Of course the Milky Way ISRF is more complicated than this monochromatic idealization, being composed of starlight photons from an array of different types of stars, of lower energy photons reprocessed by warm dust and gas and also the CMB. Energies typical for starlight photons in our galaxy are ≈ 1 − 5 eV, while photons from hot matter and the CMB are less energetic, ≈ 0.1 eV. In this case the shape of CRE spectra at a given energy reflects the combination of energy losses off of various components of the ISRF (as well as of the bulk magnetic field) and the issue of whether IC losses are dominantly Thomson-like or dominantly KN-like essentially becomes a question of the relative intensity of 0.1 eV photons to 1.0 eV photons in our local few kpc. This point has been studied in [57], where the authors discuss the impact of tuning the relative amounts of 1 eV and 0.1 eV photons in the context of an analytical model. They note that by using parameter choices that could plausibly describe our astrophysical environment, one may generically expect a "spectral pileup" (a feature that ends up looking rather more like a "bump" than a "step") in the (e + + e − ) spectrum in the neighborhood of a few hundred GeV. Their work demonstrated that one could explain any small excess that may exist in the FERMI-LAT (e + + e − ) measurement around ∼ 1 TeV while noting that the relatively large and sharp excess apparent in the ATIC (e + + e − ) measurement could not be described in this way. They also demonstrated that the putative anomaly represented in the PAMELA e + /(e + + e − ) data cannot be simultaneously explained without tuning astrophysical parameters (such as the energy density of starlight and matter density in the ISM) away from their estimated values by several orders of magnitude.
Despite the fact that tuning the ISRF does not plausibly explain the anomalous positron fraction on its own, one certainly sees that the uncertainty in the description of the ISRF can have a large effect on CRE spectra. We thus expect that a tuning of these parameters can be helpful in our present aims. In varying more parameters than just γ e we will be afforded more options with which to fit the FERMI-LAT (e + + e − ) data, each of which necessitates a particular propagation of the signal and background positrons that also enters the calculation of the e + /(e + + e − ) spectrum.
To implement this customization we use the "ISRF-Factors" that are available as parameters in the GALPROP galdef input files. The ISRF is implemented in GALPROP as a multidimensional array split into three components that represent maps of optical (starlight), FIR (reprocessed) and CMB photons as generated in the work [34]. These ISRF-Factors, referred to as F op , F F IR and F CM B here, are simply global normalizations of the three components 9 . Here, we scan the sum (F op + F F IR ) to vary the overall IC loss rate and we do not vary F CM B as there is little uncertainty in the description of the CMB. We scan the normalization of the bulk magnetic field, F B , to vary the overall synchrotron loss rate. We also vary the ratio (F op /F F IR ) to study the Thomson/KN nature of IC losses. We make no serious attempt to estimate the uncertainties in these quantities but simply limit ourselves to values within an order of magnitude of the GALPROP default set. The exact values used in this scan are given in Table 3.
Finally, we note that tuning the parameters that describe CRE sources and propagation will also affect the observed gamma-ray and radio photon spectra. The measured diffuse gamma-ray spectrum has many components: the IC and bremsstrahlung scattering of the galactic CREs, the decay of hadronically-produced π 0 mesons, known (i .e., resolved) energetic point sources, an isotropic population of γ's from extragalactic processes (such as unresolved energetic point-sources) and, in practice, contamination due to detector mis-  Table 3: Parameter scan ranges in the CRE source and radiative energy loss sector. Parameters in bold type are default (i.e. as in galdef 50p 599278, except in the case of γ e where γ e = 2.42 is default in the sense that it matches the benchmark "Model 1" of [28]).
identification of charged particles. Each of these contributions comes with its own uncertainties, whether from modeling the physical process involved or from the measurements of relevant data, that make the combination of components a non-trivial exercise.
In what follows we will show the effects of tuning these parameters on the spectrum of diffuse gamma-rays at mid-latitudes 10 . We note, however, that this calculation is approximate as the modeling of the various components that make the diffuse mid-latitude gamma-ray spectrum is still very much an active area of research. In constructing our custom astrophysical models we started with a model (galdef 50p 599278) that was designed to fit the EGRET measurement [58] of the diffuse mid-latitude gammas (see e.g., [59]), which is now known to be in tension with the more recent FERMI-LAT measurements [9] [10]. In addition, the work [60] has also presented a more accurate parameterization of the gamma spectra arising from the decays of spallatively-produced π 0 's, which is not yet incorporated into the latest public version of GALPROP. An as-yet unreleased version of GALPROP is expected to include this and other data and may also be accommodated by a new conventional astrophysical model. A recent analysis, [61], discusses the importance of uncertainties in the calculation of the diffuse mid-latitude gamma-ray spectrum in the context of constraining contributions from the annihilation of SUSY DM.

Astrophysical Scan Results
Before discussing our global fits that include the SUSY DM contributions, we can get a sense of what we have already accomplished by tuning the astrophysical parameters alone. Performing all of the steps described above leads to a set of 6615 astrophysical models. From these we take a smaller subset (524 models) to be the ones that we eventually combine with each of our ≈ 69k models from the pMSSM, so that the combined analysis could be done with a reasonable amount of CPU time. Figures 5-8 show the results for each of these models in comparison to the B/C , 10 In galactic coordinates this is 0 • ≤ l < 360 • and 10 • ≤ |b| ≤ 20 • , where the galactic longitude l is measured in the plane of the galactic disk from the sun-galactic-center axis and the galactic latitude b is the angle measured in a plane perpendicular to that of the disk, containing the sun-galactic-center axis, with the sun at its vertex.
(e + +e − ) , e + /(e + +e − ) andp/p data, respectively. In each figure we include as a benchmark the astrophysical model, "Model 1" of [28]. One can see that there is far less variation in the hadronic quantities than in the leptonic ones, as leptonic propagation parameters are far less constrained. In fact, we find that all of our 524 astrophysical models provide an excellent fit to the B/C andp/p data, that is as good, if not better, than the benchmark model. Our models also afford a good description of the (e + + e − ) data in the energy range measured by the FERMI-LAT. In the case of the positron fraction, our models have a harder spectrum at high energies compared to the benchmark case and thus have the potential to provide a better fit to the data once the supersymmetric contributions are included.
In Figure 9 we show the spectra of diffuse gamma-rays at mid-latitudes that result from these models along with the spectrum calculated with the benchmark astrophysical model. We emphasize that this is an approximate calculation, containing many uncertainties, but is nonetheless suitable for comparing the results from our custom astrophysical models with those of the benchmark. Due to these uncertainties, we do not include this data in our calculation of the global χ 2 fit. Each curve includes IC, bremsstrahlung and π 0decay contributions. One can see that our custom models exhibit more variation at low energies than they do at high energies. This is again because leptonic quantities exhibit more variation than hadronic quantities: the low-energy portion of the spectrum is dominated by gammas produced in IC up-scattering of ISRF photons by galactic CREs while the highenergy portion of the spectra is dominated by secondarily-produced π 0 -decay.
As emphasized above, we propagate our pMSSM DM signal contributions using Green's functions appropriate to each set of parameters employed in the computation of the astrophysical background spectra. Figures 10-11 demonstrate how the DM signal predictions change as we scan the astrophysical parameter space. Taking one representative SUSY model (a model with mχ0 1 = 128 GeV that annihilates almost purely to τ pairs) we display both the positron and antiproton spectra from the DM halo annihilation in each of the 524 astrophysical models. We employ the DarkSUSY default NFW profile [62] to calculate all of the signals. The choice of halo profile is found to have little impact on the positron signals, but may have a sizeable impact on the antiproton signals (as commonly used profiles differ largely only within ≈ 1 kpc of the galactic center). Nevertheless, the NFW cusp profile gives results that lie in between other commonly chosen profiles and should be a reasonable standard for our purposes.
While the benchmark "Model 1" illustrated that the FERMI-LAT (e + + e − ) measurement could be accommodated simply by choosing a relatively stiff γ e , our models show that the (e + + e − ) data can also be accommodated in diffusion models for which the (backgroundonly) e + /(e + +e − ) ratio comes significantly closer to the recent (anomalously large) PAMELA measurement. This occurs due to allowing for non-default electron energy loss parameters. We illustrate this point in Figures 12-13. Here, we color-code according to the value of resulting χ 2 for the fit in the Diffusion Zone Height versus γ e plane, first taking only the subset of astrophysical models with default IC energy loss parameters (Fig. 12) and then showing the results for all 524 models (Fig. 13). Although the best fits in each panel are found by taking Figure 5: B/C spectra for all 524 astrophysical models with no SUSY contributions (grey curves). Also shown is the B/C spectrum for the benchmark model (red curve). Data shown are from the CREAM [49], HEAO-3 [50] and ATIC experiments [51].  : (e + + e − ) spectra for all 524 astrophysical models with no SUSY contributions (grey curves). Also shown is the (e + + e − ) spectrum for the benchmark model (red curve). The FERMI (e + + e − ) data [5] is also shown. Figure 8: e + /(e + + e − ) spectra for all 524 astrophysical models with no SUSY contributions (grey curves). Also shown is the e + /(e + + e − ) spectrum for the benchmark model (red curve). Although we only fit to data from PAMELA [4] we include also data from HEAT [63], AMS01 [64], and CAPRICE94 [65], for comparison. Figure 9: Mid-latitude diffuse gamma spectra for all 524 astrophysical models with no SUSY contributions (grey curves). Also shown is the spectrum for the benchmark model (red curve). We display for comparison the EGRET [58] and FERMI [9] data, as well as preliminary FERMI data presented in [10]. Figure 10: The propagated DM positron flux (scaled by E 2 ) from a representative pMSSM model for each of our 524 astrophysical models (grey curves). Also shown is the corresponding DM positron signal from this pMSSM model propagated in the benchmark "Model 1" (red curve). The pMSSM model has an LSP with mχ0 1 = 128 GeV, which annihilates almost purely to τ pairs. Figure 11: The propagated DM antiproton flux (scaled by E 2 ) from a representative pMSSM model for each of our 524 astrophysical models (grey curves). Also shown is the corresponding DM antiproton signal from this pMSSM model propagated in the benchmark "Model 1" (red curve). The pMSSM model has an LSP with mχ0 1 = 128 GeV, which annihilates almost purely to τ pairs. a value for γ e which is as low as possible, one can see that allowing non-standard IC loss parameters widens the region of the best fit significantly. It is also apparent in Figures 12-13 that smaller diffusion zones yield markedly better fits to the data than do larger diffusion zones. We emphasize that these are astrophysics-only fits; in the next section we will show that the astrophysical models that give the best fits with SUSY contributions included have a somewhat softer γ e ≈ 2.51 − 2.55.
Another important factor in determining the global best fit is N e the normalization of primary electrons. This normalization is mostly set by the need to fit the FERMI-LAT (e + +e − ) measurement, as this data set has relatively many bins with relatively small errors as compared to the PAMELA e + /(e + +e − ) data. The primary electron normalization obviously has a large impact on the resulting positron fraction curve, and hence a large impact on the portion of the fit due to the PAMELA e + /(e + +e − ) data, especially from the low-energy data bins 11 . Even while varying other parameters, as we do here, a relatively small window of N e values yield a good fit to both data sets. This is illustrated in Figure 14, where we show the global χ 2 /dof versus N e . In figure 15 we color-code according to the primary electron normalization factor N e (rather than the global χ 2 ) in the Diffusion Zone Height versus γ e plane. We observe that the N e values corresponding to good fits to both the FERMI-LAT    (e + + e − ) and the PAMELA e + /(e + + e − ) data are found in the low γ e and small diffusion zone height corner of the plane.
We again emphasize that these results include only standard astrophysical contributions and that, with SUSY DM contributions included the primary electron normalization is allowed to float along with the boost-factor in order to determine the best fit. We now turn to a discussion of the results of the combined scan.

Results From Combined SUSY/Astrophysical Scans
In this section we describe the results of a global χ 2 to the CR data set for each combination of models in our dual scan over the astrophysical and pMSSM parameter spaces. The fit includes all (≥ 10 GeV) data bins from the FERMI-LAT (e + + e − ) measurement as well as the PAMELA e + /(e + + e − ) andp/p measurements (26, 6 and 6 bins, respectively). Errors for the FERMI-LAT data set are derived from [5] by adding statistical and systematic errors in quadrature while those for the PAMELA data sets are purely statistical, as quoted by the collaboration in [4] and [6]. The contributions to individual experimental bins from the theory spectra obtained from GALPROP/DarkSUSY are determined by numerically integrating the curves over the relevant bin. We allow the values of the boost factor, B, and of the normalization of primary electrons, N e , to float in order to obtain the best fit. This procedure leaves us with 36 degrees of freedom in calculating χ 2 /dof . We first examine the effect of a supersymmetric DM contribution to the global fit. Figure 16 shows the fraction of χ 2 /dof due to the fit to PAMELA e + /(e + +e − ) data alone as a function of the χ 2 /dof value for the global fit. In this figure, the grey points correspond to fits including both supersymmetric and astrophysical contributions for the combinations of our ∼69k pMSSM models and 524 astrophysical models for which best-fits occur with boost factors B < 200. The blue points represent the results that can be attained with astrophysical contributions alone. It is clear from the figure that the addition of supersymmetric DM can improve the fit to the data over the pure astrophysical case. The χ 2 / dof fits including SUSY DM reach a minimum of χ 2 / dof = 1.55, whereas the astrophysics-only fits obtain values reaching only as low as χ 2 / dof ≈ 2.39. Figure 16 also suggests that such a reduction of χ 2 /dof typically arises from an improvement in the fit to the PAMELA positron fraction data.
We would also like to know what values of the boost factor are required in order to obtain the best fit to the data. In Figure 17 we display the best-fit boost factor versus the reduced χ 2 for all combinations of pMSSM models and CR astrophysics models that attain their best fit when B ≤ 200. The red points in this figure represent the results from combining our ∼69k pMSSM models with the benchmark astrophysical model. We observe that the combination of our ∼69k pMSSM models with the astrophysical models in our 524 custom model set yields much improved fits with much lower boost factors compared to the benchmark astrophysical case (in fact most of the ∼69k benchmark cases do not attain Figure 16: Fraction of the full χ 2 /dof coming from the PAMELA e + /(e + + e − ) data, χ 2 pos.frac. /χ 2 tot , versus the global χ 2 /dof . The astrophysical background-only fits (blue) reach a minimum χ 2 / dof ≈ 2.39 while the fits which include a SUSY DM contribution (grey) reach a minimum χ 2 / dof ≈ 1.55. These fits were performed for each combination of our ∼69k pMSSM models and 524 astrophysical background models, though all of the grey points displayed here correspond to combinations which find their best fits for B < 200. Figure 17: Best-fit boost factor versus the global χ 2 /dof for all combinations in the astrophysics and pMSSM dual scan that attain their best fit with B < 200 (grey points). Also shown are the results from some of the fits obtained by combining our pMSSM scan with the benchmark astrophysical "Model 1" (red points) discussed in the text. Note that most of the results from the fits obtained with the benchmark model are not found in the ranges displayed here. their best fits for values of χ 2 /dof and B that fall within the ranges displayed in Fig. 17). Although this figure shows only the model combinations that have a best fit with B < 200, we allowed for boosts as high as 5000 in our numerical calculations. It is amusing to note that, although the case with the overall best fit occurs for B ≈ 3000 with χ 2 / dof = 1.544, this is hardly better than the global fits that can be obtained with B < 200, where χ 2 / dof = 1.545 is obtained for several astrophysical/SUSY model combinations with boosts in the range B = 100 − 160.
We should also remark on the perplexing feature of Figure 17 that many model combinations seem to require little or no boost to attain their best fit to the data. We find that in these cases the astrophysical background models have relatively stiff electron injection indices (γ e = 2.42 − 2.45), low radiative loss rates (F op + F F IR = 0.5 − 1) and relatively KN-like losses (F op /F F IR = 4 − 10). These features combine to give an (e + + e − ) spectrum that is very stiff, relative to other astrophysical models in our set. For these models we find that the astrophysical contribution to the positron fraction is very large in the low energy bins, leaving little room to add a flux of positrons from SUSY DM annihilation without greatly exceeding the (most tightly-constraining) positron fraction data near 10 − 20 GeV. In the end, as is evident in Figure 17, these "very-low-boost" combinations attain global fits of at best χ 2 / dof ≈ 2.0, significantly worse than the best cases, which have χ 2 / dof ≈ 1.55.
An interesting result of this analysis is that the cases with the best fit to the data set with low boost factors from among all the combinations of astrophysical and pMSSM models are dominated by both a small number of particular astrophysical models and by a small number of particular pMSSM models. This results in a few parameter sets that best describe the data, giving global fits that most frequently satisfy χ 2 / dof < 1.70 with bestfit boost factors B < 200. These astrophysical models have relatively small diffusion zones (half-heights L = 2.0 − 3.0 kpc) so that, as discussed previously, N e takes a value appropriate for fitting both the (e + + e − ) and e + /(e + + e − ) data simultaneously. These models have relatively soft injection spectra, γ e = 2.51 − 2.54, which are compensated for by IC energy losses that are more KN-like than in the default case (i.e. with F op /F F IR = 4 − 10). In addition, smaller values of the magnetic field normalization F B · B(0, 0) < ∼ 1.0µG give better fits than do larger values, though this dependence is the weakest among the parameters being scanned here.
To gauge which of our ∼69k pMSSM models provide the "best" fits to the astrophysical data, we consider which models with mχ0 1 ≥ 100 GeV most frequently (over the set of 524 astrophysical models) satisfy χ 2 / dof < 1.7 with best fit boost factors B < 200. The 10 models that most frequently satisfy these criteria are highlighted in red in Figure 18. These models are displayed in the plane comparing the ratio of DM annihilation cross sections for τ pair and bb final states versus the scaled τ pair cross section. It is clear that this set of models has a significantly larger annihilation rate to τ + τ − than to bb.
These models also have relic densities of order Ωh 2 ≈ Ωh 2 | W M AP and LSP masses that are large enough ( > ∼ 100 GeV) to make a significant impact on the CR positron fraction. We find that these models realize Ωh 2 ≈ Ωh 2 | W M AP via just the right amount of co-annihilation withτ 1 with the mass splitting: (mτ − mχ0 1 ) ∼ 15 − 20 GeV. This alone suggests that these models have sizeable τ ± annihilation rates (via t-channel exchange of the lightτ ), but with the very general mass spectra available in the pMSSM. We find that these models also annihilate with high purity into τ ± . In these models, the hadronic t-channel squark exchange is usually mass suppressed. Such high purity in the τ channel allows for a good fit with lower boost to the PAMELA positron fraction data without significantly affecting the fit to the (non-anomalous) PAMELA antiproton data. The best fit boost factors for these models are in the range B ∼ 70 − 190.
We postpone a more detailed discussion of this set of 10 pMSSM models until Section 5. We now turn to a discussion of the search for DM contributions to gamma-ray observations of dwarf satellite galaxies, where we will find that annihilation dominantly to τ -pairs is an interesting and exceptional case. Figure 18: Comparison of the τ ± and bb annihilation rates for the ten pMSSM models that most frequently (over the set of 524 astrophysical models) attain the best fits with B < 200.

Gamma-ray signals from WIMP annihilation in dwarf galaxies
Though "dark," SUSY dark matter can, via its annihilation, produce a rich gamma-ray spectrum. The detection of monochromatic γ's would represent a nearly undeniable detection of DM annihilation through loop-level channels with photons in the final state. Comparable fluxes of continuum gammas can also arise from radiation off of charged final and virtual states, as well as from the products of hadronization, fragmentation and decay (e.g., π 0 's). Such continuum spectra are highly model-dependent but may also be distinguished from astrophysical background sources by using characteristic features, such as a kinematic edge at mχ0 1 or a discernable non-power-law shape that can help to determine their DM origin upon observation.
Strategies for detecting γ's from DM annihilation are varied in their aims. If intuition from numerical simulations is correct, the inner-most part of the galaxy should house an extremely over-dense, and therefore luminous, DM cusp (along with an uncertain astrophysical background). Line signals, however, are so distinct that searches can include diffuse γ's from the entire Milky Way Halo to increase their sensitivity (and possibly overcome the loop-suppression). Another possibility is to look where there is essentially no astrophysical background so that any observation of γ's may be hypothesized to have come only from DM.
Dwarf galaxies may represent our best chance at carrying out this last type of search. Dwarfs are relatively small systems that accompany the much larger Milky Way; they harbor small collections of stars, relatively small amounts of gas and have weak magnetic fields. Their existence is anticipated from the hierarchical nature of structure formation and therefore dwarfs are expected to be highly DM-dominated (with mass-to-light ratios M/L ∼ 100 − 1000(M /L ), as determined from the kinematics of their stellar populations). Dwarfs are not expected to contain any internal astrophysical sources of γ's and though the Milky Way's own diffuse γ's present a sizable foreground, this can be subtracted reliably for dwarfs at high enough galactic latitudes. This leaves the very low background levels necessary for such low-statistics DM searches.
The full gamma-ray signal originating from a DM distribution ρ( r) factorizes into a piece dependent only on the DM distribution (or what is often called the "DM halo") and a piece dependent only on the underlying particle physics model. The number of gammas per unit area per unit time is given by: where, for a generic DM distribution ρ( r), the first factor is the integral of ρ 2 ( r) along the portion of each line-of-sight l(Ω ) for which ρ( r) = 0, over the solid angle subtended by the distribution around our detector, ∆Ω. In the second factor mχ0 1 is the LSP mass, σv is the thermal-averaged total WIMP annihilation cross-section, R is the quantity that was defined in Eqn. (2) and the quantity: is the total continuum "yield curve" with terms describing hadronization yield ("secondary" gammas), final-state radiation (FSR), virtual internal bremsstrahlung (VIB) and possible monochromatic γ's for each annihilation channel 12 . The sum runs over final state annihilation channels, indexed by i, and the weights σv B i are the annihilation rates into the i th final state (i .e., τ + τ − , bb, W + W − , etc.). We integrate the total continuum yield curve from the lower threshold E th , which is chosen according to the performance of a given detector, up to E max = m χ for all experiments considered here.
Although the calculation of yield curves is straightforward (given a particular particle physics model), estimations of the halo-dependent piece are typically fraught with great uncertainty. For instance, while much effort has been devoted to numerical simulations of DM halos similar to that of our own Milky Way Halo, the uncertainties associated with quantities such as the slope of the inner-halo cusp profile allow for the predicted flux of γ's from DM annihilating in a region near the Galactic Center (GC) to vary by more than an order of magnitude (even when specializing to fairly cuspy profiles [68]). The estimation of halo integrals for dwarf galaxies is a simpler process: here one can directly observe the stellar kinematics of the dwarf system member stars. Of course this process enjoys its own difficulties, for example in estimating the dynamical state of the dwarf DM distribution (i.e., the influence of tidal forces due to the Milky Way itself) or in determining the membership of individual stars (see e.g., [69]).
Still, the gamma-ray signals of DM annihilation incur much less uncertainty than their counterparts involving electrons and positrons. For the energies of importance here, the galaxy is essentially transparent to γ's so a non-trivial propagator is not necessary to encode their attenuation, in contrast to e ± 's. The halo integral above is carried out along a full line-of-sight because, unlike e ± 's, gammas do not diffuse but rather point directly back to their source.
Perhaps the best target for DM annihilation searches in dwarf galaxies is the ultrafaint dwarf Segue 1. Discovered relatively recently in analyses of SDSS data [70], Segue 1 is known to have a mass-to-light ratio M/L ∼ 1000(M /L ), and a preliminary analysis has been carried out to estimate its DM content using the kinematics of 24 of its member stars [71]. We will use the analysis presented in [72], appropriate to the FERMI-LAT observations of Segue 1, in which the line-of-sight halo integral was estimated by marginalizing over halo distributions, including both cusped and cored types 13 . The work [72] has recently been updated and expanded in [73].
While it may eventually be possible to obtain full spectral information, the first task at hand is to identify any excess γ's coming from dwarfs, i.e., the integrated yield from Eqns. (10) and (11). We calculate the yield curves for each of our SUSY models using the DarkSUSY package, version 5.0.4 [35], and for this analysis we do not include the contributions from γ lines 14 . We estimate the γ yield for FERMI-LAT observations of Segue 1 arising from each of our SUSY models by integrating these continuum yield curves from a low energy threshold of 100 MeV (corresponding to the FERMI-LAT threshold) up to mχ0 1 (this assumes that the LAT can dependably reconstruct γ's all the way up to our highest LSP mass, ∼ 700 GeV) and by taking the central value estimate for the Segue 1 halo integral provided in [72]. In the recent work [73] this central value estimate has been lowered significantly from the value used here. We emphasize that the results presented here are approximate and an effort to more accurately calculate the dwarf signals that can be expected from our pMSSM model set is currently underway [74].
The result is shown in Figure 19, where the horizontal line reflects the estimated sensitivity [72] for the LAT to detect a point source signal at high galactic latitudes with 3σ confidence after 1 year of operation. One can see that this sensitivity is about an order of magnitude above what is necessary to detect any of the models in our set, assuming a relic abundance consistent with thermal cosmology. In Figure 20 we examine the effect of including a boost factor to see what fraction of our model set would be excluded if effective boost factors of B = 100 or B = 1000 are incorporated. Here we are not trying to motivate such boost factors by, e.g., non-thermal cosmological history or novel halo distribution effects, rather we will use them in what follows to study the SUSY model dependence of these signals.
The vertical width of the scatter in Figures 19-20 represents the extent of the SUSY model dependence in relating σv R 2 for a given model to the total yield of gammas in our pMSSM model set. Of course there is the trivial m −2 χ 0 1 dependence apparent in Eq. (10), but it is instructive to factor this out and look for further model dependence. In this aim, we show σv R 2 /m 2 χ 0 1 as a function of mχ0 1 , emphasizing models that are NOT excluded with B = 1000 in Figure 21. In this figure the non-trivial SUSY model dependence is expressed as the extent to which models with similar σv R 2 cluster along a horizontal line. Of course, for a given boost, we expect that many models cannot be excluded simply because they have lower σv R 2 than the necessary minimum (e.g., consider the improbable case of annihilations entirely into final state γ's). However, the interesting feature of this figure is that there is a collection of ∼ 100 pMSSM models that are difficult to detect relative to others with similar values of σv R 2 (i .e., those scattered above the "exclusion line" in the Figure 19: Continuum γ yields integrated above 100 MeV as a function of σv R 2 . A point is shown for each model, orange points corresponding to pMSSM models for which Ωh 2 | LSP ≥ 0.1. The red line represents the estimated sensitivity for the LAT to detect a point source signal at high galactic latitudes with 3σ confidence [72].  1 . Models that are NOT excluded by FERMI-LAT observations of Segue 1 with B = 1000 are emphasized in dark green against the full pMSSM model set. Dark points that leak above the apparent "exclusion line" feature are pMSSM models whose continuum gamma yields are relatively difficult to detect. It is easy to see what these "difficult" models have in common by taking a look at the relative values of their annihilation rates into the various SM final states: one finds that they have especially large annihilation rates into the τ + τ − final state. For the set of models whose integrated yields are more than a factor of 2 from being excluded by FERMI observations of Segue 1 with a boost factor B = 1000, we find that in all cases the LSP annihilates into τ + τ − > ∼ 80% of the time, while the vast majority of models in this set have LSPs which annihilate into τ + τ − > ∼ 94% of the time. We highlight this set of models in Figure 22, where we show the ratio of annihilation rates into τ + τ − and bb final states as a function of the annihilation rate into the τ + τ − final state. One can see that they lie in the "leptophilic" region discussed in Section 2. Of course the focus on exclusion with B = 1000 is arbitrary and we note that an analysis of exclusion with B = 100 tells a similar story; models that are more difficult to detect have large annihilation rates to τ + τ − and, if highlighted on Figure  ( This may seem contrary to the common intuition (see e.g., [75] [7]) that large rates into the τ + τ − final state are easier to detect in γ's because of the relatively stiff π 0 decay spectrum. In Figure 23 we display the continuum γ flux spectra from the models in our set that annihilate into the τ + τ − final state > ∼ 99% of the time (red) and from those that annihilate into the final state bb > ∼ 99% of the time (blue). These models have a wide range of values for mχ0 1 and σv R 2 so we factor out this model dependence by dividing each curve by the appropriate value of σv R 2 and by using the variable x = E γ /mχ0 1 . What remains  Fig. 4. We note that these models lie in the "leptophilic" region described in the text and that the models highlighted here are similar to those highlighted in Figure 18 but with lower overall σv R 2 .
is highly model-dependent near the endpoint x → 1 where VIB can become dominant [66], but is relatively model-independent otherwise. One can see that generic τ + τ − spectra lead to higher continuum γ fluxes than do generic bb spectra, but only above a universal value x ∼ 0.15. Below this value the softer bb spectrum dominates and we see that the visibility becomes a function of the threshold energy appropriate to a given detector. The spectral crossover occurs at E γ ∼ 15 GeV for mχ0 1 = 100 GeV and at E γ ∼ 150 GeV for mχ0 1 = 1 TeV. For Cerenkov telescopes, whose threshold energies are typically E th ∼ 100 GeV, this means that the differential flux of continuum γ's from annihilation into τ + τ − would be larger than that from annihilation into the bb final state over the entire energy range above threshold. Thus, for Cerenkov telescopes, we expect that DM which annihilates dominantly into τ + τ − can be more easily detected. However, in the case of the FERMI-LAT, whose threshold energy can be taken as E th ∼ 100 MeV, the situation is quite different. The FERMI-LAT would be able to count γ's with energies well below the universal cross-over point (for mχ0 1 = 100 GeV this corresponds to integrating over the entire range shown in Figure  23), so that DM which annihilates mostly to τ + τ − would actually be more difficult for the FERMI-LAT to detect.
Several things can be said by looking deeper into the detailed mass spectra of this leptophilic model set. Starting with the identity of the nLSP, according to frequency we find:τ 1 ,ν τ ,ẽ R ,χ ± 1 andχ 0 2 nLSPs (the vast majority being of the first variety). As discussed above, it is clear that such large leptonic annihilation rates can be arranged for in models that have relatively light sleptons and relatively heavy squarks, thus suppressing or enhancing Figure 23: γ continua (calculated with DarkSUSY 5.0.4) from models in our pMSSM set for which the LSP annihilates into the τ + τ − final state greater than 99% of the time (red), and those for which the LSP annihilates into the bb final state greater than 99% of the time (blue). The curves have been rescaled according to σv R 2 and are displayed as a function of x = E γ /mχ0 1 . Note that there is a universal cross-over in the differential flux spectra near x ∼ 0.15. the associated t-channel annihilation graphs accordingly. Here, theτ 1 ,ν τ ,ẽ R nLSP cases attain leptophilia simply because theτ 1 is always relatively light (either the nLSP or close to the nLSP). In the cases where theχ ± 1 or theχ 0 2 are the nLSP, one generally expects a sizeable rate for annihilation to gauge bosons. For the model set discussed here, however, we have either that the LSPs are too light to annihilate into W or Z bosons, or that they are sufficiently purely bino-like that their couplings to gauge bosons are suppressed. These models are leptophilic simply because their colored superpartners are quite heavy (taking, e.g., the ratio of sbottom and stau masses, one finds that mb 1 /mτ 1 ∼ 2 − 7). With a relatively light (often nLSP)τ 1 one might expect that the relic densities of these models are significantly determined by stau co-annihilation (recall that models in our pMSSM set have much more generic spectra than the mSUGRA models that are usually referred to as stau co-annihilation models). One would then expect that the leptophilic models that have the largest relic density have some preferred range of mass splitting (mτ 1 − mχ0 1 ), and indeed this is the case, typically (mτ 1 − mχ0 1 ) ∼ 10 − 20 GeV. Given that the flux of cosmic-ray positrons from DM annihilation scales as the relic density squared, we expect that this specific subset of leptophilic models may be the most likely candidates for producing a sizeable excess in the positron fraction. Indeed this is the case, as we have described in Section 3.3.

Description of Best-Fit pMSSM Models
In Section 3.3 we presented the fits to the CR data fits that result from our combined scan over pMSSM and astrophysical parameter spaces. We found that the pMSSM models that most frequently provide the best fits to the data have significant annihilation rates into τ -pairs and are apparently leptophilic, as they were found to provide a significant supersymmetric DM contribution to the e + /(e + + e − ) spectrum without significant excesses over the measuredp/p data. We identified a set of ten "best" pMSSM models, which have mχ0 1 > 100 GeV and most frequently (over the set of 524 astrophysical models) satisfy χ 2 / dof < 1.7 with best fit boost factors B < 200. We now discuss these pMSSM models and their CR spectra in more detail.
We begin by presenting Table 4 where, for each of these ten models, we display the quantities that are most important for the resulting DM annihilation phenomenology. We observe that all models have LSPs with mχ0 1 ∼ 100 − 130 GeV and account for a significant fraction of the total DM density measured by WMAP. We also see that, after a thermal relic rescaling and inclusion of the appropriate boost factor, B σv R 2 ≈ 10 −24 cm 3 s −1 for all ten models. While all models in this set have significant annihilation rates into τ pairs, we note that seven of the ten annihilate into τ pairs > 80% of the time. Figures 24-26 show the e + /(e + + e − ) , (e + + e − ) andp/p spectra, respectively, for a typical case (Model 10 in Table 4) of these best-fit combinations of pMSSM and astrophysical  Table 4: Parameters describing the DM annihilation phenomenology of the 10 SUSY models which provide the best fit to CR data. For each model we display the LSP mass, mχ0 1 , the fraction of DM attributed to the LSP, via the scaling factor R (see Eqn. 2), the boost factor, B, that provides the best fit to the CR data, the boosted and scaled total annihilation cross-section B σv R 2 in units of 10 −24 cm 3 s −1 and also the ratios of annihilation rates into various final states (τ -pairs, bb, Z 0 Z 0 and W + W − ) to the total annihilation cross-section σv . models. The DM contribution is seen to produce a sizeable excess in the positron flux fraction, though still lying significantly below the trend of the PAMELA data points and appearing to fit the smaller excesses measured by HEAT and AMS-01 quite well. In the (e + + e − ) spectrum, one can see a slight bump due to the DM contribution near 100 GeV 15 . In both the (e + + e − ) andp/p spectra, we see that the pMSSM DM contribution is barely visible above the astrophysical background fluxes and provides a good fit to the data.
It is also interesting to consider the DM contributions to gamma-ray observations that arise from these pMSSM models. The relation between DM contributions to CRs and DM contributions to gamma-rays is non-trivial, however, because it is not clear how to relate the boost factor appropriate for gamma-ray quantities to that which is required to best fit the CR data. We recall that locally-measured CRE quantities give information about DM in the local ∼ kpc neighborhood of the Milky Way DM Halo, while gamma-ray observables, such as the diffuse mid-latitude gamma-ray spectrum, are composed of line-of-sight integrals of ρ 2 ( r) that extend throughout the halo. Even though the mid-latitude region is chosen in order to remove the large backgrounds along the galactic plane and to emphasize a more local diffuse gamma population, one still expects a different boost factor to apply to the gamma-ray signal (see e.g., [76]). Here, for simplicity, we will assume that the same boost  Table 4). We display curves with (black) and without (grey) the DM signal contribution. The DM contribution is seen to produce a sizeable excess. Data shown are from the PAMELA [4], HEAT [63], AMS01 [64], and CAPRICE94 [65] collaborations.  Table 4). We display curves with (black) and without (grey) the DM signal contribution. The FERMI (e + + e − ) data [5] is also shown.  Table  4). We display curves with (black) and without (grey) the DM signal contribution, though the two are indistinguishable in the figure. The PAMELAp/p data [6] is also shown.
factor that provides a best-fit the CR data should be applied to the DM contribution to the diffuse mid-latitude gamma-ray spectrum. Figure 27 shows the diffuse mid-latitude gamma-ray spectrum that can be expected from the same example (Model 10) that has been discussed in the previous figures 24-26. We calculate the DM contribution to the diffuse mid-latitude gamma-ray spectrum using the DarkSUSY default NFW halo profile. The influence of this choice on the resulting signal can be described by the dimensionless quantity (this is essentially just the first factor in Eq. 10),J where, with r = 8.5 kpc, ρ = 0.3 GeV cm −3 and ∆Ω = 2.12 for the mid-latitude region (0 • ≤ l < 360 • and 10 • ≤ |b| ≤ 20 • ), we findJ NFW = 3.21. For comparison, we find that the smooth isothermal profile implemented in DarkSUSY gives a very similar value,J ISO = 3.16, in this region.
We observe that the DM annihilation contribution to the diffuse mid-latitude astrophysical background results in a bump in the spectrum, giving a slight excess above the preliminary FERMI data for energies E ∼ 20 − 50 GeV. We emphasize again that we use the same boost as that which best fits the CR data in calculating the DM contribution to this quantity and that this is a non-trivial assumption. These diffuse mid-latitude data sets were NOT included in the simultaneous fit of the CR data in our combined scan of pMSSM and astrophysical parameters (due to the uncertainties inherent in this calculation as discussed at the end of Section 3.2.2). For the small number of pMSSM models considered here, we have repeated the global fit incorporating the diffuse mid-latitude data by including the highest-energy data bins in the preliminary FERMI data set. We find little change in the results, as the CR data sets dominate the determination of the best-fit boost factor and the resulting value of χ 2 / dof. We point out that if the appropriate boost for this quantity were a factor of 2 lower than that for the CR data, then the resulting excess due to DM would be consistent with the preliminary data at the 1σ level for all high-energy bins.  Table 4). We display curves with (black) and without (grey-solid) the DM signal contribution, as well as the DM spectrum in isolation (grey-dashed). The DM signal is calculated using the DarkSUSY default NFW halo profile. Data shown are published data from EGRET [58] and FERMI [9], as well as the preliminary FERMI data presented in [10].
It is also interesting to include the boost factors that best fit the CR data sets in computing the gamma-ray yields from DM annihilation in dwarf galaxies. We do this in the context of the FERMI-LAT observations of Segue 1 by simply multiplying the yields for each of the models described in Table 4 by the corresponding boost factors. This is again a non-trivial assumption, and in principle a very different one that has been made in the diffuse mid-latitude gamma case; the systematics involved in modeling the DM distribution of a dwarf are different than those involved in modeling the DM distribution of the Milky   [72]. We also display the ratio of the annihilation rate into τ -pairs over the total annihilation rate.
Way halo on much larger scales. The results are displayed in Table 5. We observe that the expected DM annihilation signals range from a factor of ∼ 5 above the estimated sensitivity to a factor of ∼ 2 below, and that, in keeping with the results of Section 4, the models which annihilate most purely to τ -pairs are the most difficult to detect 16 .
The result that the DM models which provide the best fits to the CR data sets should give gamma-ray contributions that are comparable to the current experimental sensitivities, to within a factor of a few, provides an exciting opportunity to confirm or refute these contributions to the anomalous CR data with measurements and analyses that will be carried out in the very near term.
We also investigate the direct detection signatures that can be expected from the pMSSM models in Table 4. Here, again, there is some ambiguity in determining the boost factor that should be applied to the DM direct detection signals, although the fact that these signals arise from DM scattering instead of DM annihilation necessitates further discussion.
While we have remained agnostic as to the origin of the boost factor throughout our analysis, simply defining it as the factor by which σv R 2 must deviate from the canonical value that would result from thermal cosmological evolution in order to best fit the CR data, the origin of the boost factor has a significant impact on the resulting direct detection rates. For example, in scenarios where the annihilation signal is enhanced due to the presence of a resonance [14], there will be no corresponding boost to the direct detection signal, which arises via t-channel scattering diagrams. In scenarios which employ non-thermal cosmological evolution [18], there will be an excess relic density of DM over that which would be expected from the usual thermal relation between σv and Ωh 2 | LSP , and the direct detection signal should be boosted accordingly. If one imagines the boost factor to originate entirely from non-trivial substructure in the local DM distribution, then scattering in direct detection experiments may or may not incur a boost, as the overdensities responsible for enhanced annihilation may or may not coincide with terrestrial detectors.
In light of this, we present Figures 28-29, where we highlight the direct detection cross sections for the ten models in Table 4 against the rest of our pMSSM model set, without any boost factors included. We note that the cross sections predicted for these models span several orders of magnitude, with some being on the edge of discovery. Figure 28: Proton spin-independent elastic scattering cross sections. Un-boosted cross sections for the 10 models described in Table 4 are highlighted (red points) against the cross sections from the entire ∼ 69k pMSSM model set. Neutron spin-independent scattering cross sections are similar and are not presented here.
Lastly, in Figures 30a-f and 31a-d we display the full superparticle mass spectrum for each of the ten pMSSM models, described in Table 4, that provide the best-fit to the CR data with B ≤ 200. The most important feature common to all models in this set is the lightτ 1 ; this being necessary for annihilations dominantly to τ -pairs. All LSPs in this set are (mostly-bino) mixtures of bino and higgsino gauge eigenstates but, we observe that the LSPs which are most purely bino-like undergo annihilations most purely into the τ + τ − final state. Figure 29: Proton spin-dependent elastic scattering cross sections. Un-boosted cross sections for the 10 models described in Table 4 are highlighted (red points) against the cross sections from the entire ∼ 69k pMSSM model set. Neutron spin-dependent scattering cross sections are similar and are not presented here.

Discussion and Conclusions
In this paper we have investigated the ability of recently measured astrophysical data to either constrain or shed light on the nature of particle dark matter in supersymmetric models. To this end we have employed the large set of phenomenologically viable MSSM models generated in [19]. We have used the software package DarkSUSY [35] to calculate the products of DM annihilations using the full details (masses, mixing angles, etc.) that describe each pMSSM model. We have stressed the need to consider the uncertainty inherent in determining astrophysical backgrounds to DM signals, and have presented an analysis in which we used the software package GALPROP [31][32] to account for these uncertainties by scanning over phenomenologically viable CR astrophysics models. We performed a simultaneous fit to the PAMELA e + /(e + + e − ) andp/p data sets as well as to the FERMI (e + + e − ) data for all possible combinations of pMSSM models and astrophysical background models, using both DarkSUSY and GALPROP to consistently combine signal and background CR fluxes.
These models are also unique in the sense that they are particularly difficult to discover or constrain from measurements of the total flux of gamma-rays coming from quiet astrophysical objects, such as the ultra-faint dwarf Segue 1 used as an example here. The contributions of these models to the diffuse mid-latitude gamma-ray spectrum yields a small excess over the preliminary FERMI data in the high energy bins, but is most likely within the uncertainties inherent in this calculation. In addition, we computed the DM direct detection rates for these pMSSM models and found that they span several orders of magnitude, with some models being on the edge of discovery. The CR fluxes generated by this class of pMSSM models stand out because they are produced by DM annihilations which are dominantly into the τ + τ − final state, a property which is due to having very light (often nLSP)τ 1 and relatively heavy colored superpartner mass spectra. Models in this class which give rise to the largest annihilation fluxes, i .e. the largest σv R 2 , have mass splittings (mτ 1 − mχ0 1 ) ∼ 10 − 20 GeV, which are not so small that they excessively dilute the relic density, R, and not so large that they excessively decrease the total cross-section, σv .
Finally we remark on a few caveats concerning our results and on prospects for future observations or constraints on DM.
We first note that while we produced simultaneous fits to the e + /(e + + e − ),p/p and (e + + e − ) data with sizeable excesses over background in e + /(e + + e − ) and no corresponding excess inp/p, the quality of the best fits is not as good as is perhaps desired. It is clear that the shape of the DM contribution to the e + /(e + + e − ) fraction from this class of "best" pMSSM models is not a perfect fit to the sharp rise observed by PAMELA, although it would easily explain the more modest excesses observed in the AMS-01 and HEAT data. There is also the question of how to obtain a boost factor in the range B ∼ 70 − 190, which is a number that is probably still too large to be explained by the current uncertainty in the normalization or substructure of the DM halo distribution of the Milky Way [77] [13]. Other works have attempted the explanation of such factors in terms of non-standard cosmology, or non-trivial modifications to the BSM particle model, while here we have remained strictly agnostic as to the origin of the boost, seeking only to quantify its size.
We also note that more prosaic astrophysical explanations (e.g. a local population of pulsars [78]) could account for the rise in the positron fraction spectrum observed by PAMELA. In order to confirm a SUSY DM explanation from scenarios such as those discussed in this paper, one would likely need to see that subsequent measurements of the positron fraction suggest an electroweak scale (∼ 100 − 200 GeV) edge feature, and that the upward slope of the positron fraction is not as steep as that observed by PAMELA. Such assertions would require an experiment with exceptional proton/positron ID, even for particles of energy in the ∼ 100 − 200 GeV range, and large enough effective volume to beat down statistical errors in such higher-energy bins. There is at least one near-future experiment, AMS-02, that may be capable of providing such a measurement on the ∼ 5 year time scale 17 .
In the same time frame the LHC is expected to provide its first exploration of the TeV scale. If the BSM particle physics describing Nature is similar to the class of models discussed above, with mχ0 1 ∼ 100 − 200 GeV, a (typically nLSP)τ 1 of mass mτ 1 ≈ mχ0 1 + (10 − 20) GeV, and relatively heavy colored superpartners, one may expect to observe many events with energetic τ leptons and / E T . Conversely, observations of such events at the LHC may be used to inform indirect detection strategies or, with non-observation of DM signals in γ and CR fluxes, to provide an indirect bound on the relic density of the observedχ 0 1 . Such a limit may be an important complement to the constraints that may be attained via direct detection, especially in the case of a leptophilic LSP.
One may also expect that the astrophysical uncertainties discussed here may gradually decrease as data from experiments such as FERMI [80], AMS-02 [81], CREAM [82] and TRACER [83] are expected to enable a host of analyses that may improve and refine CR propagation models, though one certainly expects the case of CR e ± , whose transport depends strongly on the detailed features of the local galactic neighborhood, to remain quite a challenge.