CosRayMC: a global fitting method in studying the properties of the new sources of cosmic e ± excesses
aa r X i v : . [ a s t r o - ph . C O ] F e b CosRayMC: a global fitting method in studying the properties ofthe new sources of cosmic e ± excesses Jie Liu a , Qiang Yuan b , Xiao-Jun Bi b , Hong Li a,c , and Xinmin Zhang a,ca Theoretical Physics Division, Institute of High Energy Physics,Chinese Academy of Science, P.O.Box 918-4, Beijing 100049, P.R.China b Key Laboratory of Particle Astrophysics,Institute of High Energy Physics, Chinese Academy of Science,P.O.Box 918-3, Beijing 100049, P.R.China and c Theoretical Physics Center for Science Facilities (TPCSF),Chinese Academy of Science, Beijing 100049, P.R.China
Recently PAMELA collaboration published the cosmic nuclei and electron spectrawith high precision, together with the cosmic antiproton data updated, and theFermi-LAT collaboration also updated the measurement of the total e + e − spectrumto lower energies. In this paper we develop a Markov Chain Monte Carlo (MCMC)package CosRayMC , based on the GALPROP cosmic ray propagation model to studythe implications of these new data. It is found that if only the background electronsand secondary positrons are considered, the fit is very bad with χ ≈ .
68. Takinginto account the extra e + e − sources of pulsars or dark matter annihilation we cangive much better fit to these data, with the minimum χ ≈ .
83. This means theextra sources are necessary with a very high significance in order to fit the data.However, the data show little difference between pulsar and dark matter scenarios.Both the background and extra source parameters are well constrained with thisMCMC method. Including the antiproton data, we further constrain the branchingratio of dark matter annihilation into quarks B q < .
5% at 2 σ confidence level. Thepossible systematical uncertainties of the present study are discussed. PACS numbers: 95.35.+d,96.50.S-
I. INTRODUCTION
A very interesting progress made in the recent years in cosmic ray (CR) physics is thediscovery of the excesses of positrons and electrons by several space- and ground-based ex-periments [1–5]. The positron and electron excesses challenge the traditional understandingof CR background. Many theoretical models were proposed to explain these new phenom-ena, including the astrophysical scenarios (e.g., [6–11], and see [12] for a review) and thedark matter (DM) scenario (e.g., [13–16])Most recently Fermi-LAT team reported the updated measurements of the total spec-trum of electrons and positrons, extending to energies as low as several GeV [17]. PAMELAcollaboration updated the observation of the ¯ p/p ratio and the absolute antiproton flux [18],and reported the measurement of the pure electron spectrum at the first time [19]. Theantiproton data extend to 180 GeV, without any hint of deviation from the backgroundcontribution [18]. For the PAMELA electron data, although there is no significant spectralfeature above 30 GeV other than a single power-law, it is also consistent with models includ-ing extra e + e − sources to explain the positron excess [19]. There is also hint of hardening ofthe electron spectrum compared with the low energy part ( <
20 GeV), although the solarmodulation may be important at these low energies.Since the accumulation of high-quality CR data, it is now important to extract moreinformation from these data, i.e., estimating the CR background and the possible extrasource parameters. Previously, one always constrained one or two parameters with otherparameters fixed. This may lead to biased results, especially, when the parameters arestrongly correlated. The global fitting procedure searches the maximum likelihood in themultiple dimensional parameter space rather than a reduced one and can give the posteriordistribution by marginalizing other parameters in the Bayesian approach. The Markov ChainMonte Carlo (MCMC) procedure, whose computational time scales approximately linearlywith the number of parameters, makes it possible to survey in a very large parameter spacewith the least computational cost.In our previous studies [20, 21], we employed the MCMC method to fit the parameters ofthe DM scenario as well as the background parameters proposed to explain the e + e − excesses.However, the propagation of CRs is treated with a semi-analytical way following Refs. [22–24]. A more precise description of the CR propagation is given by the numerical models,such as GALPROP [25] and DRAGON [26], in which most of the relevant physical processesare taken into account, and the realistic astrophysical inputs like the interstellar medium(ISM) and the interstellar radiation field (ISRF) are adopted. There are many parameters inthe CR propagation model, and it is very difficult to have a full and systematical survey ofthe parameter space. After embedding the numerical CR propagation tool into the MCMCsampler, we can use it to constrain the model parameters in a more efficient way. Recentlythere are also several other works using the MCMC method to study the CR propagation[27–29].In this work, we develop a package CosRayMC (Cosmic Ray MCMC), which is comprisedof GALPROP, PYTHIA [30] and MCMC sampler, to revisit the models to explain thepositron and electron excesses and derive the constraints on the model parameters with thelatest data. Two kinds of the extra e + e − sources are considered: the pulsar scenario andthe DM annihilation scenario. In both scenarios we consider the continuous distributionof the sources, although it is possible that one or several nearby pulsars or DM subhalosmay explain the excesses [6, 7, 31–33]. As pointed out in [34], it was unlikely that the fluxfrom any single pulsar was significantly larger than that from others given a large number ofknown nearby and energetic pulsars. For DM subhalos the location and mass are uncertain,which also makes the constraints of the model parameters difficult.This paper is organized as follows. We briefly introduce the propagation of Galactic CRsin Sec. II. The description of the methodology is given in Sec. III. The fitting results ofthe pulsar and DM scenarios are presented in Sec. IV. Finally we give the conclusion anddiscussion in Sec. V. II. PROPAGATION OF GALACTIC COSMIC RAYS
The charged particles propagate diffusively in the Galaxy due to the scattering withrandom magnetic field. There are interactions between the CR particles and the ISM and/orthe ISRF, which will lead to fragmentation, catastrophic or continuous energy losses of theparticles. For unstable nuclei the radioactive decay also needs to be taken into account.In addition, the overall convection driven by the stellar wind and reacceleration due tothe interstellar shock will also affect the distribution function of CRs. For each species ofparticles we have a partial differential equation to describe the propagation process, withthe general form [35] ∂ψ∂t = Q ( x , p )+ ∇· ( D xx ∇ ψ − V c ψ )+ ∂∂p p D pp ∂∂p p ψ − ∂∂p h ˙ pψ − p ∇ · V c ψ ) i − ψτ f − ψτ r , (1)where ψ is the density of cosmic ray particles per unit momentum interval, Q ( x , p ) is thesource term, D xx is the spatial diffusion coefficient, V c is the convection velocity, D pp isthe diffusion coefficient in momentum space used to describe the reacceleration process, ˙ p ≡ d p/ d t is the momentum loss rate, τ f and τ r are time scales for fragmentation and radioactivedecay respectively. Solving the partially coupled equations for all kinds of particles, we canget the propagated results of the CR spectra and spatial distributions. For more detailsabout the terms in Eq. (1) please refer to the recent review paper [35].The secondary-to-primary ratios such as B/C and (Sc+Ti+V)/Fe, and the unstable-to-stable ratios of secondary particles such as Be/ Be and Al/ Al, are often used toconstrain the propagation parameters because the ratios can effectively avoid the sourceparameters. Then one can use the spectra of the primary particles to derive the sourceparameters. There are some studies to constrain the propagation parameters based on thecurrently available data (e.g., [25, 36, 37]). In [27–29] the MCMC method was adopted tofit both the propagation and source parameters of CRs. However, due to the quality of theobservational data, the constraints on the propagation parameters are not very effective,and there may be also large systematic uncertainties.Since we focus on the electron and positron data in this work, which are not very effectiveto constrain the propagation parameters, we adopt the best fitting results of the propagationparameters given in [29]. The main propagation parameters are: D = 6 . × cm s − , δ = 0 . v A = 39 . − , z h = 3 . γ n = 1 .
91 (below 10 GV) and γ n = 2 .
40 (above 10 GV). What we should keep in mind isthat there might be systematic errors of the determination of the propagation parameters(see e.g., [38]), and the quantitative results of this work might be affected.We note that the newest measurements of the CR proton and Helium spectra by CREAM[39] and PAMELA [40] showed remarkable deviation of single power-law spectra. Further-more the spectral indices of proton and Helium are different. Such new features challenge thetraditional understanding of the CR origin, acceleration and propagation (see e.g. [41, 42]for some explanations). In the present work we keep the study in the traditional frame anddo not try to reproduce the detailed structures of the data. We find that the above adoptedinjection parameters with proper adjusting of the normalization can basically reproduce thePAMELA-CREAM proton data, as shown in the left panel of Fig. 1. However, when com-paring with the Helium data, the model expectation seems to be systematically lower. Weincrease the relative abundance of Helium in GALPROP by 30%, and give roughly consis-tent result with PAMELA data (right panel of Fig. 1). The CREAM data can still not bereproduced. Since the contribution to the secondary particles from CR Helium is only asmall fraction, we do not think a higher Helium flux within a factor of 2 will substantiallyaffect our following study. -1 E k F l u x ( G e V / nu c l eon m - s - s r - ) E k (GeV/nucleon) φ = 500 MVLIS galpropAMSBESSATIC2PAMELACREAM 10 -1 -1 E k F l u x ( G e V / nu c l eon m - s - s r - ) E k (GeV/nucleon) φ = 500 MVLIS galpropBESSCAPRICE94ATIC2PAMELACREAM FIG. 1: Proton (left) and Helium (right) spectra of the GALPROP model calculation, comparedwith the observational data. In each panel the higher line labelled “LIS” represents the localinterstellar flux, and the lower one is the flux after the solar modulation. References of the dataare: proton — AMS [43], BESS [44], ATIC2 [45], PAMELA [40], CREAM [39]; Helium — BESS[44], CAPRICE94 [46], ATIC2 [45], PAMELA [40], CREAM [39].
III. METHODOLOGYA. CosRayMC
The
CosRayMC (Cosmic Ray MCMC) code is built up by embedding the CR propagationcode GALPROP (GALPROP v51) into the MCMC sampling scheme. The MCMC techniqueis widely applied to give multi-dimensional parameter constraints from observational data.Following Bayes’ Theorem, the posterior probability of a model (which we refer to a seriesof parameters ~θ ) given the data (which are described by D ) is P ( ~θ | D ) ∝ L ( ~θ ) P ( ~θ ) , (2)where L ( ~θ ) = P ( D | ~θ ) is the likelihood function of the model ~θ for the data D , and P ( ~θ ) is theprior probability of the model parameters. In performing CosRayMC , MCMC is employedto generate a random sample from the posterior distribution P ( ~θ | D ) which are fair samplesof the likelihood surface. Based on the sample, we can get the estimate of the mean values,the variance as well as the confidence level of the model parameters. For details please referto [47].The key improvement, compared with our previous study, is that the calculation of thelikelihood L ( ~θ ) are given by calling GALPROP to simulate the mock observations. By doingso, we can have a more precise description of the CR propagation. And as the better data-set is provided in the future, our code can also be used to make a determination of the CRpropagation parameters.For the part with DM contribution to the CRs, the PYTHIA (PYTHIA v6.4) simula-tion code is employed to calculate the final spectra of electrons, positrons and antiprotons[30]. Such spectra are then injected into the Galaxy and propagated with GALPROP. ThePYTHIA code is also embedded in the CosRayMC code.
B. Parameters
We assume the injection spectrum of the background electrons to be a broken power-lawfunction with spectral indices γ / γ below/above E br . Note that for the shock accelerationscenario, the injection spectrum of particles can not be too hard [48]. We set the priorsthat γ > . γ > . A bkg , taken as the flux of electrons at 25 GeV, is also regarded as a free parameter.For the background positrons and antiprotons, we adopt the GALPROP model predictedresults with the best fitting source and propagation parameters given in [29]. Consideringthe fact that there are uncertainties about the ISM density distribution, the hadronic inter-action model, and the propagation parameters determined from the secondary-to-primaryratio data, we will further employ factors c e + and c ¯ p to rescale the absolute fluxes of thesesecondary particle.For energies below ∼
30 GeV the solar modulation effect is important and needs tobe considered. In this work we adopt the force-field approximation to calculate the solarmodulation [49]. The modulation potential depends on the solar activity. For the periodwhich PAMELA works the modulation potential is estimated to be φ = 450 −
550 MV [40].In our MCMC fit, the modulation potential φ is also taken as a free parameter. Thus, forthe background model we have 7 parameters in total P bkg = { γ , γ , E br , A bkg , φ, c e + , c ¯ p } . (3)Pulsars are thought to be the most natural candidates to generate high energy positronsand electrons through the cascade of electrons accelerated in the magnetosphere [31, 50].The spectrum of e + e − escaped from the pulsars can be parameterized as power-law with acutoff at E c , d N/ d E ∝ A psr E − α exp( − E/E c ), where the power-law index α ranges from 1to 2.2 according to the radio and gamma-ray observations [31]. The cutoff energy of theinjected e + e − ranges from several tens GeV to higher than TeV, depending on the modelsand parameters of the pulsars [34, 50]. For the spatial distribution of pulsars, we adopt thefollowing form [51] f ( R, z ) ∝ (cid:18) RR ⊙ (cid:19) a exp (cid:20) − b ( R − R ⊙ ) R ⊙ (cid:21) exp (cid:18) − | z | z s (cid:19) , (4)where R ⊙ = 8 . z s ≈ . a = 2 .
35 and b = 5 .
56. There will be a furthernormalization factor A psr .Alternatively, DM annihilation or decay models are widely employed to explain the e + e − excesses. Since to date one can not distinguish DM annihilation from decay with the e + e − data [20, 21], we only consider the annihilation scenario here. Considering the non-excessof PAMELA antiproton data, the DM annihilation final states need to be lepton-dominated[15, 16]. Similar as the previous works [20, 21], we do not employ the annihilation finalstates based on any DM models, but to constrain the final states from the data we use amodel-independent way instead. Such constraints will be helpful for the understanding ofthe DM particle nature.The annihilation final states are assumed to be two-body e + e − , µ + µ − , τ + τ − and q ¯ q , with Since the antiproton productions of various quark flavors do not differ much from each other, here we donot distinguish quark flavors but adopt u ¯ u channel as a typical one of quark final states. branching ratios B e , B µ , B τ and B q respectively. It was also shown that the interactionswith low-mass intermediate bosons which then decay to lepton pairs could fit the CR e + e − data well [52, 53]. However, since the spectral shapes of the electrons and positrons fromthese models do not have distinct properties which can be easily distinguished from µ + µ − and τ + τ − , we do not include such final states in this study. The resulting positron, electronand antiproton spectra of these final states are calculated using the PYTHIA simulationpackage [30].The DM density profile is assumed to be an Einasto type [54] ρ ( r ) = ρ − exp (cid:20) − α (cid:18) r α r α − − (cid:19)(cid:21) , (5)where α ≈ . r − ≈ . ρ − ≈ .
14 GeV cm − according to the recent highresolution simulation Aquarius [55]. The corresponding local density of DM is about 0 . − , which is consistent with the results of a larger local density derived in recentstudies [56–58]. Then the source function of positron produced from DM annihilation is q ( E, r ) = h σv i m χ X f B f d N d E (cid:12)(cid:12)(cid:12)(cid:12) f × ρ ( r ) , (6)where m χ is the mass of the DM particle, h σv i is the velocity weighted annihilation crosssection, d N d E (cid:12)(cid:12) f is the positron production rate from channel f of one annihilation.The most general parameter space is summerized as below. P tol = {P bkg } , background {P bkg , A psr , α, E c } , pulsar {P bkg , c ¯ p , m χ , h σv i , B e , B µ , B τ , B u } , DM . (7)For the DM annihilation scenario, we further consider two cases. One is to use the e + e − data only in the fit, which can be directly compared with the pulsar scenario. In the othercase, the antiproton data are included in order to constrain the coupling to quarks of DM.In the case without the antiproton data, B u and c ¯ p are set to zero. IV. RESULTSA. Background
First we run a fit for the pure background contribution. The data included in the fitare PAMELA positron fraction [1], PAMELA electron [19], Fermi-LAT total e + e − [17], andHESS total e + e − [3, 4].For the PAMELA positron fraction data, we only employ the data points with E > e + e − motion[61] can give consistent description to the PAMELA and previous data. The demodulatedinterstellar spectra of electrons and positrons are also consistent with the conventional CRbackground model expectation. Here we do not consider the detailed solar modulationmodels. But note that the solar modulation model may affect the quantitative fitting results.This is one kind of systematical errors.The best-fit results of the positron fraction and electron (or e + e − ) spectra are shownin Fig. 2. The best-fit parameters are compiled in Table I, with the minimum χ ≈ . /
113 = 3 .
68. Such a large reduced χ means that the fit is far from acceptable butsystematics dominated. Using a χ goodness-of-fit test we find that the background onlycase is rejected with a very high significance ∼ σ . That is to say the data strongly favorthe existence of additional degrees of freedom. From Fig. 2 we can clearly see that thebackground model under-estimates both the positron fraction and the total e + e − fluxes.The extra sources of e + e − are necessary to explain simultaneously the positron fraction andtotal e + e − data.The poor fit makes nonsense to discuss the physical implication of the parameters. There-fore we leave the discussion of the parameters in the next subsections, when the extra sourcesof e + e − are taken into account.0 -2 -1 e + / ( e - + e + ) E (GeV)bkg10 -2 -1 e + / ( e - + e + ) E (GeV) AMSHEAT94+95HEAT00PAMELA 10 E d N / d E ( G e V m - s - s r - ) E (GeV)bkg e - bkg e + +e - E d N / d E ( G e V m - s - s r - ) E (GeV) PAMELAATICHESSFermi-LAT
FIG. 2: Best-fit positron fraction (left) and electron spectra (right) of the pure background case.Note that in the right panel the PAMELA data are for pure electrons, while other data are forthe sum of electrons and positrons. References of the observational data are: positron fraction —AMS [62], HEAT94+95 [63], HEAT00 [64], PAMELA [1]; electron — PAMELA [19], ATIC [2],HESS [3, 4], Fermi-LAT [17].
B. Pulsar scenario
Compared with the conventional GALPROP background model, γ = 1 . ± .
08 isconsistent with 1 .
60 as given in [65]. Note that there is a correlation between parameters γ and the modulation potential φ , as shown in Fig. 3. The high energy injection spectrum γ = 2 . ± .
02 is softer than 2 .
54 as adopted in the conventional model [65]. In theconventional model only the background contribution is employed to fit the data, whilehere an additional component from the extra sources are added together. Thus the highenergy spectrum can be much softer. This might be important for the understanding of thebackground contribution to the CR electrons.There is a factor c e + ≈ . pp interaction to generate positrons. Compared with the pp collision modelBadhwar77 [67] as adopted in GALPROP, the Kamae06 model gives systematically fewerpositrons, especially for energies from several to tens of GeV [8]. If these input uncertainties1 TABLE I: Fitting parameters with 1 σ uncertainties or 2 σ limits. Note that for the “bkg” casethe reduced χ is too large that the uncertainties of the parameters should not be statisticallymeaningful. bkg bkg+pulsar bkg+DM (without ¯ p ) bkg+DM (with ¯ p ) γ < . C.L. ) 1 . +0 . − . . +0 . − . < . C.L. ) γ . ± .
008 2 . ± .
021 2 . ± .
018 2 . ± . A bkg ) − . ± . − . ± . − . ± . − . ± . E br ( GeV ) 3 . +0 . − . . +0 . − . . +0 . − . . +0 . − . φ ( GV ) 0 . ± .
019 0 . ± .
071 0 . +0 . − . . ± . c e + . ± .
037 1 . ± .
111 1 . +0 . − . . ± . c ¯ p — — — 1 . ± . A psr ) — − . +0 . − . — — α — 1 . +0 . − . — — E c — 0 . +0 . − . — — M χ ( T eV ) — — 2 . +0 . − . . +0 . − . log[ σv ( cm s − )] — — − . ± . − . ± . B e — — < . C.L. ) < . C.L. ) B µ — — < . C.L. ) < . C.L. ) B τ — — 0 . +0 . − . . +0 . − . B u — — — < . C.L. ) χ /d.o.f .
676 0 .
833 0 .
827 1 . A bkg is in unit of cm − sr − s − M eV − A psr is in unit of cm − sr − s − M eV − are better understood in the future, we can include the background positron productionwith a more general form in the MCMC fit.The fitted parameters for pulsars are α = 1 . +0 . − . , E c = 0 . + . − . TeV. The con-straints of the pulsar parameters are weaker when compared with the background param-eters. This is because the cutoff energy is mainly determined by HESS data. However,for HESS data, the very large systematic errors due to poor absolute energy calibration2 P r obab ili t y P r obab ili t y P r obab ili t y P r obab ili t y P r obab ili t y P r obab ili t y E c γ E b r φ α γ E c γ br φ α FIG. 3: Two dimensional constraints of part of the model parameters for the pulsar scenario.In off-diagonal figures, the inner contour denotes the 68% confidence level (C.L.) while the outercontour stands for 95% C.L. The diagonal figures are one dimensional probability distributions forcorresponding parameters. (not shown in Fig. 4) makes it impossible to precisely determine E c . In the calculationthe relative systematic uncertainty of the flux is estimated to be (Γ − E/E [17]. For ∼
15% uncertainty of the energy scale, the flux uncertainty is ∼
30% and ∼
45% for energiesbelow and above ∼ E c and α , A psr , the other twoparameters are also relatively less constrained.3 -2 -1 e + / ( e - + e + ) E (GeV)bkgpulsartotal10 -2 -1 e + / ( e - + e + ) E (GeV) AMSHEAT94+95HEAT00PAMELA 10 E d N / d E ( G e V m - s - s r - ) E (GeV)bkg e - bkg e + +e - pulsartotal10 E d N / d E ( G e V m - s - s r - ) E (GeV) PAMELAATICHESSFermi-LAT
FIG. 4: Same as Fig. 4 but for the model with pulsars as the extra sources of e + e − . C. Dark matter annihilation scenario
First we consider the case without including the antiproton data. We found that thebackground parameters are similar with the pulsar scenario. The overall χ value of thebest-fit is about 90 .
2, which is very close to χ ≈ .
6. That is to say both the astrophysicalscenario and DM scenario can give comparable fit to the e + e − data. We may not be ableto discriminate these two scenarios from the e + e − spectra [34, 68–71], but need to resort toother probes such as the e + e − anisotropy [72] and γ -rays [73].The fitted mass and cross section of DM are m χ = 2 . +0 . − . TeV, h σv i = 4 . +1 . − . × − cm s − . Similar as the pulsar case, these parameters are also less constrained due tothe large uncertainties of the HESS data. We can also determine the branching ratios todifferent flavors of leptons: B e (2 σ ) < . B µ (2 σ ) < .
38 and B τ = 0 . +0 . − . . These resultsare qualitatively consistent with our previous fit in [20], but with larger uncertainties becausewe include the systematic errors of HESS data in the present work. We should be cautiousthat the determination of the branching ratios may suffer from systematic uncertaintiesof the background electrons and positrons. As discussed in the previous sub-section, thepresent understanding of the background positrons is far from precise. Therefore the quotedfitting uncertainties of the branching ratios may be under-estimated.Then we consider the case with antiproton data. The coupling to quarks of DM an-nihilation is taken into account. The two dimensional constraints of part of parametersare present in Fig. 5. Compared with the case without antiproton data, we find that thebackground parameters are slightly different. The reason for the change is mainly due to4 P r obab ili t y P r obab ili t y P r obab ili t y P r obab ili t y −29.4 P r obab ili t y log[ σ v/2M χ ] γ E b r ( G e V ) φ ( G V ) γ l og [ σ v / M χ ] γ br (GeV)3.544.5−29.4−29.2 φ (GV)0.30.40.5−29.4−29.2 FIG. 5: Same as Fig. 3 but for the DM scenario with antiproton data. the low energy spectrum of the background antiprotons. As shown in Fig. 6 the calculatedsecondary (including tertiary) antiprotons are basically consistent with the PAMELA datafor energies higher than several GeV, but under-estimate the flux for the low energy part.This problem was pointed out years ago [77, 78], that in the diffusive reacceleration modelwhich gives the best fit to the B/C data the antiprotons are under-estimated. There seemsto be a contradiction between the B/C and the antiproton data . To better understand Note, however, in [79] the calculated antiproton flux based on the propagation parameters fitted accordingto the B/C data [36] was consistent with the observational data. As pointed out in [78] the fit in [36] -6 -5 -4 -3 -2 -1 -1 F l u x ( G e V - m - s - s r - ) E k (GeV)10 -6 -5 -4 -3 -2 -1 -1 F l u x ( G e V - m - s - s r - ) E k (GeV) AMSBESS00BESS02PAMELA FIG. 6: Fitted 2 σ range of the antiproton fluxes of the DM annihilation scenario. References ofthe data are: AMS [74], BESS00 [75], BESS02 [76], PAMELA [18]. this issue we may need more precise measurement about the B/C data. Back to this work,a lower antiproton flux at low energies will require a smaller solar modulation potential,and therefore the background parameters γ , γ and E br will change accordingly due to thecorrelations among them as shown in Fig. 5.One thing important is that we can derive a self-consistent constraint on the quarkbranching ratio of DM annihilation. The marginalized 2 σ upper limit of B u is about 0 . B u is several times smaller. Thisis probably due to the bad fit to the antiproton data with the background model. We haverun a test to employ the background antiproton spectrum used in [21] and found that the 2 σ upper limit of B u is about 2 . B u is systematics dominated. Better understanding of thebackground contribution to the antiproton flux is necessary to further address this issue.We plot the 2 σ range of the antiproton fluxes in Fig. 6. was actually based on the high energy data. And in their results the low energy B/C was in fact overestimated. Furthermore the antiproton production cross sections would also lead to uncertainties [78].Alternatively, in [80] a unified model to explain the B/C and antiproton data was proposed, with anempirical modification of the diffusion coefficient at low speed of particles. V. CONCLUSION AND DISCUSSION
Recently more and more observational data of CRs with unprecedented precision areavailable, which makes it possible to better approach the understanding of the basic problemsof CRs. Based on MCMC analysis in our previous work [20, 21], we embed GALPROP andPYTHIA into MCMC sampler and study the implication of the newest CR data, includingthe positron fraction, electrons (pure e − and e + e − ) and antiprotons from PAMELA, Fermi-LAT and HESS experiments in this work.We work in the frame of diffusive reacceleration propagation model of CRs. The prop-agation parameters are adopted according to the fit to currently available B/C data [29],with a slight adjustment of the Helium abundance to better match the PAMELA data [40].We find that the pure background to explain the CR e + e − data is disfavored with a veryhigh significance. Therefore it is strongly implied that we may need some extra sources toproduce the positrons/electrons.We then consider two different scenarios, which are widely discussed in recent literature,to explain the e + e − excesses. One is the astrophysical scenario with pulsars as the typicalexample, and the other is the DM annihilation scenario. Using the global fitting method,we can determine the parameters of both the background and the extra sources. The lowenergy spectral index of the background electrons is consistent with the conventional modeladopted in the previous study, while the high energy index is softer in this work. We findthat, together with the background contribution, both of these scenarios can give very goodfit to the data. For the case with only the e + e − data, the goodness-of-fit of the pulsarscenario and DM scenario are almost identical. With the antiproton data included, we canfurther constrain the coupling to quarks of the DM model. The branching ratio to quarksof the DM annihilation final states is constrained to be < .
5% at 2 σ confidence level.The constraint on the quark branching ratio (or the cross section to quarks) allows usto make a comparison with the direct detection experiments, assuming some effective in-teraction operators between DM and standard model (SM) particles. As shown in [81],generally the constraint of DM-SM coupling from direct experiments is much stronger thanthe indirect search for spin-independent interactions. However, for the spin-dependent in-teractions, the indirect search constraint could be comparable or stronger than the directexperiments. In Fig. 7 we give the constraints on the spin-dependent DM-neutron scattering7cross section according to the cross section of DM annihilation to quarks, assuming the axialvector interaction form of DM particles and SM fermions. The ratio between h σv i and σ SD χn is 5 c (1 + m n /m χ ) /m n P q q − m q /m χ m q [81], where c is light speed, m n is the mass ofneutron and m q is the mass of quark. The sum is for all the flavors of quarks. Consideringthe uncertainties of the background antiprotons, we show the results of two antiproton back-ground models with shaded region: the result calculated in this work and the one used in[21]. The results show that the constraint on the spin-dependent cross section between DMand nucleon from PAMELA antiproton data is much stronger than the direct experiments. −41 −40 −39 −38 −37 −36 m χ (GeV) σ S D χ n ( c m ) PAMELA pbarCDMS−2006XENON10−2008
FIG. 7: Constraints on the spin-dependent DM-neutron scattering cross section from the PAMELAantiproton data assuming axial vector interaction form of DM particles and SM fermions. Theshaded region represents the uncertainties of the astrophysical background of antiprotons. Alsoshown are the results from CDMS [82] and XENON10 [83].
Although the existence of the extra sources is evident, we still can not discriminate thepulsar and DM scenarios right now. These two models both can fit the data very well, with χ = 0 .
833 and χ = 0 .
827 respectively. Other new probes such as the e + e − anisotropyand γ -rays are needed to distinguish these models.The γ -rays should give further constraints on the models, especially for the DM annihi-lation scenario. However, the diffuse γ -ray data of Fermi-LAT are still in processing, andthe modeling of the γ -ray data seems not trivial. Furthermore the γ -rays are more sensitive8to the density profile of DM, which is not the most relevant parameter of the present study.Therefore we do not include the γ -ray data. It should be a future direction to include the γ -ray data in the CosRayMC package.There are some systematical uncertainties of the current study, which are mainly due tothe lack of knowledge about the related issues and can be improved in the future. First,the solar modulation which affects the low energy spectra of CR particles ( E .
30 GeV)if not clear. In this work we use the simple force-field approximation [49] to deal withthe solar modulation. But such a model seems to fail to explain the low energy data ofthe PAMELA positron fraction, which may indicate the charge dependent modulation effect[59–61]. Second, the propagation model is adopted as the diffusive reacceleration model withthe parameters best fitting the B/C data [29]. There are uncertainties of the propagationparameters. Furthermore the diffusive reacceleration model may also have some systematicalerrors when comparing with all of the CR data. For example the antiprotons below severalGeV is under-estimated in this model. We need to understand the propagation model betterwith more precise data. Third, the secondary positron and antiproton production may sufferfrom the uncertainties from the ISM distribution and the hadronic cross sections. All ofthese uncertainties may affect the quantitative results of the current study. Nevertheless,with more and more precise data available in the near future (e.g., AMS02, which waslaunched recently), and better control of the above mentioned systematical errors, it is aright direction to do the global fit to derive the constraints and implication on the models.
Acknowledgements
We thank Hong-Bo Hu and Zhao-Huan Yu for helpful discussion. The calculation istaken on Deepcomp7000 of Supercomputing Center, Computer Network Information Centerof Chinese Academy of Sciences. This work is supported in part by National Natural ScienceFoundation of China under Grant Nos. 90303004, 10533010, 10675136, 10803001, 11033005and 11075169, by 973 Program under Grant No. 2010CB83300, by the Chinese Academy ofScience under Grant No. KJCX2-EW-W01 and by the Youth Foundation of the Institute9of High Energy Physics under Grant No. H95461N. [1] O. Adriani, G. C. Barbarino, G. A. Bazilevskaya, R. Bellotti, M. Boezio, E. A. Bogomolov,L. Bonechi, M. Bongi, V. Bonvicini, S. Bottai, et al., Nature , 607 (2009), 0810.4995.[2] J. Chang, J. H. Adams, H. S. Ahn, G. L. Bashindzhagyan, M. Christl, O. Ganel, T. G. Guzik,J. Isbert, K. C. Kim, E. N. Kuznetsov, et al., Nature , 362 (2008).[3] F. Aharonian, A. G. Akhperjanian, U. Barres de Almeida, A. R. Bazer-Bachi, Y. Becherini,B. Behera, W. Benbow, K. Bernl¨ohr, C. Boisson, A. Bochow, et al., Phys. Rev. Lett. ,261104 (2008), 0811.3894.[4] F. Aharonian, A. G. Akhperjanian, G. Anton, U. Barres de Almeida, A. R. Bazer-Bachi,Y. Becherini, B. Behera, K. Bernl¨ohr, A. Bochow, C. Boisson, et al., Astron. Astrophys. ,561 (2009).[5] A. A. Abdo, M. Ackermann, M. Ajello, W. B. Atwood, M. Axelsson, L. Baldini, J. Bal-let, G. Barbiellini, D. Bastieri, M. Battelino, et al., Phys. Rev. Lett. , 181101 (2009),0905.0025.[6] H. Y¨uksel, M. D. Kistler, and T. Stanev, Phys. Rev. Lett. , 051101 (2009), 0810.2784.[7] D. Hooper, P. Blasi, and P. Dario Serpico, Journal of Cosmology and Astro-Particle Physics , 25 (2009), 0810.1527.[8] T. Delahaye, R. Lineros, F. Donato, N. Fornengo, J. Lavalle, P. Salati, and R. Taillet, Astron.Astrophys. , 821 (2009), 0809.5268.[9] H.-B. Hu, Q. Yuan, B. Wang, C. Fan, J.-L. Zhang, and X.-J. Bi, Astrophys. J. Lett. ,L170 (2009), 0901.1520.[10] N. J. Shaviv, E. Nakar, and T. Piran, Physical Review Letters , 111302 (2009), 0902.0376.[11] P. Blasi, Physical Review Letters , 051104 (2009), 0903.2794.[12] Y.-Z. Fan, B. Zhang, and J. Chang, International Journal of Modern Physics D , 2011(2010), 1008.4646.[13] L. Bergstr¨om, T. Bringmann, and J. Edsj¨o, Phys. Rev. D , 103520 (2008), 0808.3725.[14] V. Barger, W.-Y. Keung, D. Marfatia, and G. Shaughnessy, Phys. Lett. B , 141 (2009),0809.0162.[15] M. Cirelli, M. Kadastik, M. Raidal, and A. Strumia, Nuclear Physics B , 1 (2009), , 023512 (2009), 0811.0176.[17] M. Ackermann, M. Ajello, W. B. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri,B. M. Baughman, K. Bechtol, F. Bellardi, et al., Phys. Rev. D , 092004 (2010).[18] O. Adriani, G. C. Barbarino, G. A. Bazilevskaya, R. Bellotti, M. Boezio, E. A. Bogomolov,L. Bonechi, M. Bongi, V. Bonvicini, S. Borisov, et al., Physical Review Letters , 121101(2010), 1007.0821.[19] O. Adriani, G. C. Barbarino, G. A. Bazilevskaya, R. Bellotti, M. Boezio, E. A. Bogomolov,M. Bongi, V. Bonvicini, S. Borisov, S. Bottai, et al., Physical Review Letters , 201101(2011).[20] J. Liu, Q. Yuan, X. Bi, H. Li, and X. Zhang, Phys. Rev. D , 023516 (2010), 0906.3858.[21] J. Liu, Q. Yuan, X. J. Bi, H. Li, and X. M. Zhang, ArXiv e-prints:0911.1002 (2009), 0911.1002.[22] J. Lavalle, Q. Yuan, D. Maurin, and X. Bi, Astron. Astrophys. , 427 (2008), 0709.3634.[23] J. Lavalle, J. Pochon, P. Salati, and R. Taillet, Astron. Astrophys. , 827 (2007).[24] D. Maurin, R. Taillet, and C. Combet, ArXiv Astrophysics e-prints (2006), astro-ph/0609522.[25] A. W. Strong and I. V. Moskalenko, Astrophys. J. , 212 (1998), arXiv:astro-ph/9807150.[26] C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, Journal of Cosmology and Astro-ParticlePhysics , 18 (2008), 0807.4730.[27] A. Putze, L. Derome, D. Maurin, L. Perotto, and R. Taillet, Astron. Astrophys. , 991(2009), 0808.2437.[28] A. Putze, L. Derome, and D. Maurin, Astron. Astrophys. , A66+ (2010), 1001.0551.[29] R. Trotta, G. J´ohannesson, I. V. Moskalenko, T. A. Porter, R. Ruiz de Austri, and A. W.Strong, Astrophys. J. , 106 (2011), 1011.0037.[30] T. Sj¨ostrand, S. Mrenna, and P. Skands, Journal of High Energy Physics , 26 (2006),arXiv:hep-ph/0603175.[31] S. Profumo, ArXiv e-prints: 0812.4457 (2008), 0812.4457.[32] D. Hooper, A. Stebbins, and K. M. Zurek, Phys. Rev. D , 103513 (2009), 0812.3202.[33] M. Kuhlen and D. Malyshev, Phys. Rev. D , 123517 (2009), 0904.3378.[34] D. Malyshev, I. Cholis, and J. Gelfand, Phys. Rev. D , 063005 (2009), 0903.1310.[35] A. W. Strong, I. V. Moskalenko, and V. S. Ptuskin, Annual Review of Nuclear and Particle Science , 285 (2007), arXiv:astro-ph/0701517.[36] D. Maurin, F. Donato, R. Taillet, and P. Salati, Astrophys. J. , 585 (2001).[37] M. Pato, D. Hooper, and M. Simet, J. Cosmol. Astropart. Phys. , 22 (2010), 1002.3341.[38] D. Maurin, A. Putze, and L. Derome, Astron. Astrophys. , A67+ (2010), 1001.0553.[39] H. S. Ahn, P. Allison, M. G. Bagliesi, J. J. Beatty, G. Bigongiari, J. T. Childers, N. B. Conklin,S. Coutu, M. A. DuVernois, O. Ganel, et al., Astrophys. J. Lett. , L89 (2010), 1004.1123.[40] O. Adriani, G. C. Barbarino, G. A. Bazilevskaya, R. Bellotti, M. Boezio, E. A. Bogomolov,L. Bonechi, M. Bongi, V. Bonvicini, S. Borisov, et al., Science , 69 (2011), 1103.4055.[41] Y. Ohira and K. Ioka, Astrophys. J. Lett. , L13+ (2011), 1011.4405.[42] Q. Yuan, B. Zhang, and X.-J. Bi, ArXiv e-prints (2011), 1104.3357.[43] J. Alcaraz, B. Alpat, G. Ambrosi, H. Anderhub, L. Ao, A. Arefiev, P. Azzarello, E. Babucci,L. Baldini, M. Basile, et al., Phys. Lett. B , 27 (2000).[44] T. Sanuki, M. Motoki, H. Matsumoto, E. S. Seo, J. Z. Wang, K. Abe, K. Anraku, Y. Asaoka,M. Fujikawa, M. Imori, et al., Astrophys. J. , 1135 (2000), arXiv:astro-ph/0002481.[45] A. D. Panov, J. H. Adams, Jr., H. S. Ahn, K. E. Batkov, G. L. Bashindzhagyan, J. W. Watts,J. P. Wefel, J. Wu, O. Ganel, T. G. Guzik, et al., Bulletin of the Russian Academy of Science,Phys. , 494 (2007), arXiv:astro-ph/0612377.[46] M. Boezio, P. Carlson, T. Francke, N. Weber, M. Suffert, M. Hof, W. Menn, M. Simon, S. A.Stephens, R. Bellotti, et al., Astrophys. J. , 457 (1999).[47] A. Lewis and S. Bridle, Phys. Rev. D , 103511 (2002), arXiv:astro-ph/0205436.[48] M. A. Malkov and L. O’C Drury, Reports of Progress in Physics , 429 (2001).[49] L. J. Gleeson and W. I. Axford, Astrophys. J. , 1011 (1968).[50] L. Zhang and K. S. Cheng, Astron. Astrophys. , 1063 (2001).[51] D. R. Lorimer, in Young Neutron Stars and Their Environments , edited by F. Camilo &B. M. Gaensler (2004), vol. 218 of
IAU Symposium , p. 105.[52] L. Bergstr¨om, J. Edsj¨o, and G. Zaharijas, Physical Review Letters , 031103 (2009),0905.0333.[53] I. Cholis, G. Dobler, D. P. Finkbeiner, L. Goodenough, and N. Weiner, Phys. Rev. D ,123518 (2009), 0811.3641.[54] J. Einasto, Trudy Inst. Astrofiz. Alma-Ata , 87 (1965).[55] J. F. Navarro, A. Ludlow, V. Springel, J. Wang, M. Vogelsberger, S. D. M. White, A. Jenkins, C. S. Frenk, and A. Helmi, Mon. Not. Roy. Astron. Soc. , 21 (2010), 0810.1522.[56] R. Catena and P. Ullio, J. Cosmol. Astropart. Phys. , 4 (2010), 0907.0018.[57] M. Pato, O. Agertz, G. Bertone, B. Moore, and R. Teyssier, Phys. Rev. D , 023531 (2010),1006.1322.[58] P. Salucci, F. Nesti, G. Gentile, and C. Frigerio Martins, Astron. Astrophys. , A83+(2010), 1003.3101.[59] J. M. Clem, D. P. Clements, J. Esposito, P. Evenson, D. Huber, J. L’Heureux, P. Meyer, andC. Constantin, Astrophys. J. , 507 (1996).[60] B. Beischer, P. von Doetinchem, H. Gast, T. Kirn, and S. Schael, New Journal of Physics ,105021 (2009).[61] P. Bobik, M. J. Boschini, C. Consolandi, S. Della Torre, M. Gervasi, D. Grandi, K. Kudela,S. Pensotti, and P. G. Rancoita, ArXiv e-prints (2010), 1011.4843.[62] M. Aguilar, J. Alcaraz, J. Allaby, B. Alpat, G. Ambrosi, H. Anderhub, L. Ao, A. Arefiev,P. Azzarello, L. Baldini, et al., Phys. Lett. B , 145 (2007), arXiv:astro-ph/0703154.[63] S. W. Barwick, J. J. Beatty, A. Bhattacharyya, C. R. Bower, C. J. Chaput, S. Coutu, G. A.de Nolfo, J. Knapp, D. M. Lowder, S. McKee, et al., Astrophys. J. Lett. , L191 (1997),arXiv:astro-ph/9703192.[64] S. Coutu, A. S. Beach, J. J. Beatty, A. Bhattacharyya, C. Bower, M. A. Duvernois,A. Labrador, S. P. McKee, S. Minnick, D. Muller, et al., in International Cosmic Ray Con-ference (2001), vol. 5 of
International Cosmic Ray Conference , p. 1687.[65] A. W. Strong, I. V. Moskalenko, and O. Reimer, Astrophys. J. , 962 (2004), arXiv:astro-ph/0406254.[66] T. Kamae, N. Karlsson, T. Mizuno, T. Abe, and T. Koi, Astrophys. J. , 692 (2006),arXiv:astro-ph/0605581.[67] G. D. Badhwar, R. L. Golden, and S. A. Stephens, Phys. Rev. D , 820 (1977).[68] V. Barger, Y. Gao, W. Y. Keung, D. Marfatia, and G. Shaughnessy, Phys. Lett. B , 283(2009), 0904.2001.[69] D. Chowdhury, C. J. Jog, and S. K. Vempati, ArXiv e-prints (2009), 0909.1182.[70] M. Pohl, Phys. Rev. D , 041301 (2009), 0812.1174.[71] M. Pato, M. Lattanzi, and G. Bertone, J. Cosmol. Astropart. Phys. , 20 (2010), 1010.5236.[72] I. Cernuda, Astroparticle Physics , 59 (2010), 0905.1653. [73] J. Zhang, X. J. Bi, J. Liu, S. M. Liu, P. F. Yin, Q. Yuan, and S. H. Zhu, Phys. Rev. D ,023007 (2009), 0812.0522.[74] M. Aguilar, J. Alcaraz, J. Allaby, B. Alpat, G. Ambrosi, H. Anderhub, L. Ao, A. Arefiev,P. Azzarello, E. Babucci, et al., Phys. Rept. , 331 (2002).[75] Y. Asaoka, Y. Shikaze, K. Abe, K. Anraku, M. Fujikawa, H. Fuke, S. Haino, M. Imori,K. Izumi, T. Maeno, et al., Physical Review Letters , 051101 (2002), arXiv:astro-ph/0109007.[76] S. Haino, K. Abe, H. Fuke, T. Maeno, Y. Makida, H. Matsumoto, J. W. Mitchell, A. A.Moiseev, J. Nishimura, M. Nozaki, et al., in International Cosmic Ray Conference (2005),vol. 3 of
International Cosmic Ray Conference , pp. 13–+.[77] I. V. Moskalenko, A. W. Strong, J. F. Ormes, and M. S. Potgieter, Astrophys. J. , 280(2002).[78] I. V. Moskalenko, A. W. Strong, S. G. Mashnik, and J. F. Ormes, Astrophys. J. , 1050(2003), arXiv:astro-ph/0210480.[79] F. Donato, D. Maurin, P. Salati, A. Barrau, G. Boudoul, and R. Taillet, Astrophys. J. ,172 (2001).[80] G. di Bernardo, C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, Astroparticle Physics ,274 (2010), 0909.4548.[81] J.-M. Zheng, Z.-H. Yu, J.-W. Shao, X.-J. Bi, Z. Li, and H.-H. Zhang, ArXiv e-prints (2010),1012.2022.[82] D. S. Akerib, M. S. Armel-Funkhouser, M. J. Attisha, C. N. Bailey, L. Baudis, D. A. Bauer,P. L. Brink, P. P. Brusov, R. Bunker, B. Cabrera, et al., Phys. Rev. D , 011102 (2006),arXiv:astro-ph/0509269.[83] J. Angle, E. Aprile, F. Arneodo, L. Baudis, A. Bernstein, A. Bolozdynya, L. C. C. Coelho,C. E. Dahl, L. Deviveiros, A. D. Ferella, et al., Physical Review Letters101