Multimessenger constraints on the dark matter interpretation of the Fermi-LAT Galactic center excess
MMultimessenger constraints on the dark matter interpretation of the
Fermi -LATGalactic center excess
Mattia Di Mauro
Istituto Nazionale di Fisica Nucleare, via P. Giuria, 1, 10125 Torino, Italy ∗ Martin Wolfgang Winkler
Stockholm University and The Oskar Klein Centre for Cosmoparticle Physics, Alba Nova, 10691 Stockholm, Sweden † An excess of γ rays in the data measured by the Fermi
Large Area Telescope in the direction ofthe Galactic center has been reported in several publications. This excess, labeled as the Galacticcenter excess (GCE), is detected analyzing the data with different interstellar emission models,point source catalogs and analysis techniques. The characteristics of the GCE, recently measuredwith unprecedented precision, are all compatible with dark matter particles (DM) annihilating inthe main halo of our Galaxy, even if other interpretations are still not excluded. We investigatethe DM candidates that fit the observed GCE spectrum and spatial morphology. We assume asimple scenario with DM annihilating into a single channel but we inspect also more complicatedmodels with two and three channels. We perform a search for a γ -ray flux from a list of 48 MilkyWay dwarf spheroidal galaxies (dSphs) using state-of-the-art estimation of the DM density in theseobjects. Since we do not find any significant signal from the dSphs, we put upper limits on theannihilation cross section that result to be compatible with the DM candidate that fits the GCE.However, we find that the GCE DM signal is excluded by the AMS-02 ¯ p flux data for all hadronicand semi-hadronic annihilation channels unless the vertical size of the diffusion halo is smaller than2 kpc – which is in tension with radioactive cosmic ray fluxes and radio data. Furthermore, AMS-02 e + data rule out pure or mixed channels with a component of e + e − . The only DM candidatethat fits the GCE spectrum and is compatible with constraints obtained with the combined dSphsanalysis and the AMS-02 ¯ p and e + data annihilates purely into µ + µ − , has a mass of 60 GeV androughly a thermal cross section. PACS numbers:
I. INTRODUCTION
Several groups have discovered an excess in the γ -raydata collected by the Fermi
Large Area Telescope (
Fermi -LAT) in the direction of the Galactic center (see, e.g.,[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]). This sig-nal, called the Galactic center excess (GCE), has beendetected using different background models, constitutedby the flux of point and extended sources, interstellaremission,
Fermi bubbles and an isotropic component, andby performing the analysis with different data selectionand analysis techniques. The GCE has a spectral en-ergy distribution (SED, measured as E dN/dE in unitsof GeV/cm /s/sr) that peaks at a few GeV, a spatialmorphology that is roughly spherically symmetric and itscentroid is located in the dynamical center of the MilkyWay [9, 11, 14].The origin of the GCE is still a mystery. Refs. [15, 16],by applying wavelet analysis and non-Poissonian tem-plate fitting techniques to Fermi -LAT data, derived com-pelling evidence for the existence of a faint populationof sources located in the Galactic center with proper-ties that can explain the GCE. The presence of these ∗ Electronic address: [email protected] † Electronic address: [email protected] sources could be interpreted as a population of millisec-ond pulsars located around the Galactic bulge. Theseresults are supported by Refs. [17, 18] that modeled theGalactic stellar bulge, a possible tracer of pulsars in thecenter of the Galaxy, using a nuclear bulge and a boxybulge template. They demonstrate that fitting the GCEwith these two templates they obtain a much better fitthan using a DM model. This result implies that theGCE is not spherically symmetric since the model usedin Refs. [17, 18] has a boxy shape.Very recently, Refs. [19, 20] have shown that the non-Poissonian template fitting method can misattribute un-modeled point sources or imperfections in the IEM to asignal of a faint population of sources or DM. These re-sults cast serious doubts on the robustness of the resultspresented in [16] and the conclusion that the GCE is dueto a population of pulsars. In addition, Ref. [21] has ap-plied wavelet analysis, similarly to what has been done in[15], to about 10 years of
Fermi -LAT data using the latest4FGL catalog released by the
Fermi -LAT Collaboration[22]. They find that the GCE is still present but they donot find any compelling evidence for the existence of afaint population of un-modeled sources.Outbursts of cosmic rays (CRs) from the Galactic cen-ter have been proposed as possible interpretations for theGCE (see, e.g., [23, 24, 25]). In these alternative scenar-ios the GCE is explained by γ rays produced through in-verse Compton scattering (ICS) of high-energy electrons a r X i v : . [ a s t r o - ph . H E ] J a n and positrons on the interstellar radiation field (ISRF)photons or by CR protons interacting with the interstel-lar gas and producing π , which subsequently decays into γ rays. These mechanisms, however, provide γ -ray sig-nals not fully compatible with the GCE properties. Forexample, the hadronic scenario (i.e., CR protons) pre-dicts a γ -ray signal that is distributed along the Galacticplane, since the π decay process is correlated with thedistribution of gas present in the Milky Way disk [24].Instead, a leptonic outburst would lead to a signal thatis approximatively spherically symmetric but it requiresa complicated scenario with at least two outbursts to ex-plain the morphology and the intensity of the excess.Very recently, Ref. [14] has provided the most preciseresults for the GCE properties yet. They confirm thatthe GCE SED is peaked at a few GeV and has a high en-ergy tail significantly detected up to about 50 GeV. TheSED changes in normalization by roughly 60% when us-ing different interstellar emission models (IEMs), data se-lections and analysis techniques. The spatial distributionof the GCE is compatible with a dark matter (DM) tem-plate modeled with a generalized Navarro-Frenk-White(NFW) density profile with slope γ = 1 . − .
3. The en-ergy evolution of the GCE spatial morphology has beenstudied with unprecedented precision between 0 . − γ average value, which is 1.25, is detected. The GCEcentroid is compatible with the dynamical center of theMilky Way and its morphology is compatible with aspherical symmetric NFW profile. In particular, by fit-ting the DM spatial profile with an ellipsoid they finda major-to-minor axis ratio (aligned along the Galacticplane) between 0.8-1.2 when running the analysis withdifferent IEMs.The characteristics of the GCE published in Ref. [14]make γ rays from DM particle interactions a viable in-terpretation. In fact DM is predicted to be distributed inthe Milky Way as a spherically symmetric halo with itscentroid located in the dynamical center of the Galaxy.Moreover, the signal morphology is expected to be en-ergy independent, i.e. the value of the NFW slope ( γ )found to fit the GCE morphology data should not varywith energy. The GCE SED can be well modeled as γ rays produced by DM particles annihilating into b ¯ b witha thermal annihilation cross section [9, 10], which is theproper cross section to explain the observed density ofDM in the Universe [26]. All these characteristics makethe GCE very appealing for the DM interpretation.If DM is the origin of the GCE, γ rays should beemitted from these elusive particles also in Milky Waydwarf spheroidal galaxies (dSphs). dSphs are amongthe most promising targets for the indirect search ofDM with γ rays because gravitational observations in-dicate that they have a high DM density, i.e. a largemass-to-luminosity ratio of the order of 100 − γ rays and, as a consequence, they could provide tightconstraints on the DM interpretation of the GCE.The indirect search of DM is performed also withCR antiparticles, such as positrons ( e + ) and antiprotons(¯ p ), which are among the rarest cosmic particles in theGalaxy. e + and ¯ p fluxes have been precisely measuredby the AMS-02 experiment on the International SpaceStation up to almost 1 TeV [33, 34]. Very recently, theAMS-02 Collaboration has released the data for severalCR species with 7 years of data including new measure-ments for cosmic e + and ¯ p fluxes [35]. e + mainly originate from secondary production, dueto the spallation reactions of CRs with interstellar gasatoms, and from PWNe (see, e.g., [36, 37, 38]). No clearsignal of DM can be claimed with CR e + because of thelarge uncertainty mainly due to the possible PWN con-tribution [38]. In fact, recent observations of ICS halosdetected in γ rays around close pulsars (see, e.g., [37, 39])have provided clear evidences that PWNe inject e ± inthe interstellar space. However, it is still not clear whichfraction of pulsar spin-down energy is converted into e ± and what is exactly the acceleration process that takesplace in these sources (see, e.g., [38, 40]). The flux dataof these particles can provide very tight constraints forleptonic DM, i.e. annihilating or decaying into the e + e − , µ + µ − and τ + τ − channels, (see, e.g., [41, 42]).Different groups have found an excess of ¯ p , with re-spect to the secondary production, in the data collectedafter 4 years of mission by AMS-02 between 5-20 GeV[33]. Its significance was determined as 3 − σ dependingon the analysis technique used and the ¯ p production crosssections employed in the analysis. The excess was inter-preted in terms of DM particles with a mass of 60 − b ¯ b [43, 44, 45, 46]. Possible links tothe GCE were considered. Very recently, Ref. [47] investi-gated the presence of the ¯ p excess by fully including in theanalysis the uncertainties on the ¯ p cross sections [48], CRpropagation and correlations in the AMS-02 systematicerrors. They find that when including all these sources ofuncertainties the global significance is reduced to below1 σ . As a result the constraints for the DM interpretationof the GCE for the hadronic channels, i.e. DM annihilat-ing into quarks, might be strong also using ¯ p CR data.In this paper we investigate the DM interpretation ofthe GCE with a combined analysis of the targets that arethe most promising for the search in γ rays, i.e. the Galac-tic center and dSphs, and using the flux data of AMS-02for positrons and antiprotons which are among the rarestCRs. It is the first time ever that such an analysis for DMis performed at the same time in different astrophysicaltargets and cosmic particles and with a consistent modelfor the DM density distribution and coupling parameters.We first determine the DM density in the Galaxy using atthe same time the GCE surface brightness data reportedin [14] and the results recently obtained in Ref. [49] fromobservations of the rotation curve of the Milky Way. Wefit the GCE spectrum and find the relevant DM param-eters, mass, annihilation cross section and branching ra-tio, in case of annihilation into single, double and triplechannels. Then, we search for a DM signal in a com-bined analysis of Fermi -LAT data in the direction of 48dSphs reported in [28]. We include in the analysis theuncertainty for the DM density of these objects. Sincewe do not find any significant flux, we put upper limitsfor the annihilation cross section for the cases that bestfit the GCE spectrum. We also search for a DM signal inthe latest AMS-02 measurements of e + and ¯ p data [35].For ¯ p we use an analysis, as in Ref. [47], that accounts forthe uncertainties on the ¯ p cross sections, CR propagationand correlations between AMS-02 data. Instead, for CR e + , given the current uncertainty in the possible flux ofthese particles from PWNe, we derive constraints on aDM contribution with a conservative and an optimisticapproach. In the former we require the sum of secondarybackground and DM signal not to overshoot the AMS-02 data, while in the latter we include a possible pulsarcontribution through an analytic function similar to theapproach in Ref. [41]. Finally, we compare the DM can-didates that fit well the GCE data with the constraintsfound from dSphs, ¯ p and e + and provide the channelsand coupling parameters that satisfy all the above obser-vations.The paper is organized as follow: in Sec. II we explainthe model we use for the calculation of the γ -ray, ¯ p and e + flux from DM and for the secondary production. InSec. III we estimate the DM density in the Galaxy us-ing the GCE data and latest rotation curve data of theMilky Way. We fit the GCE with γ rays from DM usingdifferent annihilation models reporting the best-fit valuesof the relevant DM coupling parameters. In Sec. IV weperform a DM search from dSphs in the Fermi -LAT dataand produce limits for the annihilation cross section thatwe compare with the DM candidates compatible with theGCE SED. In Sec. V and VI we will investigate the com-patibility of the DM interpretation of the GCE with ¯ p and e + flux data from AMS-02. Finally, in Sec. VII wedraw our conclusions. II. MODEL FOR COSMIC PARTICLEPRODUCTION FROM DARK MATTERA. γ -ray flux from dark matter
1. Prompt emission
The γ -ray emission from DM particle interactions isusually calculated including two components. The firstone is the so-called prompt emission that is due to thedirect production of γ rays through an intermediate an-nihilation channel. The prompt emission is calculated as follows: dNdE = 12 r (cid:12) π (cid:18) ρ (cid:12) M DM (cid:19) ¯ J × (cid:104) σv (cid:105) (cid:88) f Br f (cid:18) dN γ dE (cid:19) f , (1)where M DM is the DM mass, ρ (cid:12) is the local DM density, r (cid:12) is the distance of the Earth from the center of theGalaxy, (cid:104) σv (cid:105) defines the annihilation cross section timesthe relative velocity, averaged over the Galactic velocitydistribution function. Moreover, ¯ J is the geometricalfactor averaged over the viewing solid angle ∆Ω of ourregion of interest (ROI) that is 40 ◦ × ◦ centered in theGalactic center as in Ref. [14].( dN γ /dE ) f is the γ -ray spectrum from DM annihila-tion for a specific annihilation channel labeled as f and Br f is its branching ratio. We take ( dN γ /dE ) f fromRef. [50] where this quantity has been calculated us-ing the Pythia Monte Carlo code (version 8.162). Inparticular we consider the tables reported at this web-page for the case where electroweak corrections are also in-cluded.¯ J is calculated as the integral performed along the lineof sight (l.o.s., s ) of the squared DM density distribution ρ divided for ∆Ω:¯ J = 1∆Ω (cid:90) ∆Ω d Ω (cid:90) l.o.s. dsr (cid:12) (cid:18) ρ ( r ( s, Ω)) ρ (cid:12) (cid:19) . (2)We parametrize ρ with a generalized NFW (gNFW) DMdensity function [51]: ρ gNFW = ρ s (cid:16) rr s (cid:17) γ (cid:16) rr s (cid:17) − γ , (3)or with an Einasto profile [52]: ρ Einasto = ρ s exp (cid:18) − α (cid:20)(cid:18) rr s (cid:19) α − (cid:21)(cid:19) , (4)or using a cored Burkert profile [53]: ρ Burkert = ρ s (cid:16) rr s (cid:17) (cid:20) (cid:16) rr s (cid:17) (cid:21) . (5)The parameters ρ s and r s are the normalization and scaleradius of the DM density profile, which has to be foundcalibrating ρ on the observed distribution of DM in theGalaxy. The results will be given with the gNFW andEinasto profiles since, as we will demonstrate in Sec. III Athe Burkert profile is not adequate to fit the rotationcurve and the GCE surface brightness data.
2. Inverse Compton scattering emission
In case DM particles annihilate into leptonic channels,i.e. e + e − , µ + µ − and τ + τ − , there is a secondary pro-duction of γ rays that becomes relevant. This involves e ± produced from the prompt emission and that subse-quently generate γ rays through ICS on the ISRF pho-tons. This component is particularly relevant for theGalactic center where the density of the starlight anddust components of the ISRF are roughly a factor of 10higher than the local one (see, e.g., [54]).The flux of γ rays for ICS at energy E is calculated as[50, 55, 56]: dNdE ( E ) = r (cid:12) π (cid:18) ρ (cid:12) M DM (cid:19) (cid:90) ∆Ω d Ω (cid:90) l.o.s. dsr (cid:12) ×× (cid:90) M DM E dE e N e ( E e , r ( s, Ω)) P ( E, E e , r ( s, Ω)) , (6)where N e ( E e , r ( s, Ω)) is the density of e ± produced withenergy E e from DM at a position r and P ( E, E e , r ) isthe power of γ rays produced for ICS on the ISRF. P isdefined as: P ( E, E e , r ) = 3 σ T c m e c E e (cid:90) / (4 γ ) dq (cid:18) − m e c qE e (1 − (cid:15) ) (cid:19) × n ( (cid:15) ( q, r )) G ( q ) q , (7)where σ T is the Thomson cross section, m e is the electronrest mass, n ( (cid:15) ) is the ISRF spectrum with photon energy (cid:15) , Γ = 4 (cid:15)γ/ ( m e c ) and q = (cid:15)/ (Γ( γm e c − (cid:15) )). G ( q ) iscalculated from the Klein-Nishina cross section is definedas: G ( q ) = 2 q log q + (1 + 2 q )(1 − q ) + ψ (1 − q )2(1 − ψ ) , (8)where ψ = E/E e .In order to find N e we solve the equation for the prop-agation of e ± in the Galactic diffusive halo. The prop-agation for e ± , that is dominated by energy losses forICS and synchrotron radiation and the diffusion on theirregularities of the Galactic magnetic field, is modeledas: ∂ t N e − ∇ · { K ( E ) ∇N e } + ∂ E { b ( E ) N e } = Q ( E, r ) , (9)where b ( E ) represents the energy losses, K ( E ) the dif-fusion, and Q ( E, r ) the source term for the productionof e ± from DM. Other processes usually taken into ac-count for CR nuclei are negligible for the propagation of e − (see, e.g., [57]). Assuming homogeneous energy lossesand diffusion in the Galaxy, the solution of the propaga-tion equation is found as: N e ( E e , r ) = (cid:90) M DM E e dE s (cid:90) dV d r G ( E e , r ← E s , r s ) Q ( E s , r s ) , (10)where G ( E e , r ← E s , r s ) is the Green function which ac-counts for the probability that e ± emitted at an initialGalactic position r s and with an energy E s is detected ata final position r and energy E e . Since the boundaries of the propagation zone do virtually not affect the solutionfor e ± in the Galactic center region, we can employ thefree Green function: G ( E e , r ← E s , r s ) = 1 b ( E e )( πλ ) / exp (cid:18) − ( r − r s ) λ (cid:19) . (11) λ is the propagation length for e ± affected by energylosses and diffusion: λ ( E, E s ) = 4 (cid:90) E e E s dE (cid:48) K ( E (cid:48) ) b ( E (cid:48) ) . (12)The source term Q ( E s , r s ) for DM is calculated as: Q ( E s , r s ) = (cid:18) ρ ( r s ) ρ (cid:12) (cid:19) (cid:88) f Br f (cid:18) dN e dE s (cid:19) f , (13)where ( dN e /dE s ) f is the spectrum of e ± produced fromDM particle interactions and it depends on the specificannihilation channel assumed and labeled in the equationwith f .In our ROI centered in the Galactic center, we canneglect diffusion because the propagation of e ± is dom-inated by energy losses since the starlight and infraredcomponents of the ISRF are a factor of about 30 and 8larger than in the local Galaxy respectively [54]. There-fore, for the calculation of the ICS γ -ray flux produced inour ROI we can neglect diffusion. Assuming that K ( E ) isparametrized as K ( E ) = K E δ , the typical timescale ofdiffusion is calculated as τ ∼ L /K ( E ) ∼ · E − δ Myrwith δ ≈ .
40 and K = 3 · cm /s [47]. Instead, theenergy losses in the Galactic center region have a char-acteristic time scale τ ∼ E/ ( b ( E )) ∼ · E − . Myr. At10 GeV the energy loss τ is thus a factor of about 100smaller than the one for diffusion confirming thus thatwe can neglect diffusion. In this scenario, Eq. 6 for theICS flux simplifies to the following expression [58]: dNdE = r (cid:12) π (cid:18) ρ (cid:12) M DM (cid:19) ¯ J × (cid:104) σv (cid:105) (cid:88) f Br f ×× (cid:90) M DM E e dE e P ( E, E e ) Y f ( E e ) b ( E e ) , (14)where Y f ( E e ) is defined as Y ( E e ) = (cid:82) M DM E e ( dN e /dE e ) f .In order to demonstrate further that diffusion can beneglected, we calculate the γ -ray emission for ICS in-cluding and neglecting diffusion. We assume the ISRFmodel for the Galactic center as in Ref. [54] and theparametrization of diffusion as in Ref. [47]. We performthe calculation for two leptonic channels µ + µ − and τ + τ − and the hadronic channel b ¯ b and for a DM mass and crosssection of 50 GeV and 3 × − cm /s. These are roughlythe DM parameters that best fit the GCE spectrum (seeSec. III B 1). We show in Fig. 1 the result of this calcula-tion. The case with the µ + µ − channel is the one, amongthe three shown, for which the difference between the casewith and without diffusion is more evident in the totalflux because the ICS component gives the largest contri-bution with respect to the prompt emission. Instead, forthe b ¯ b channel since the ICS component is negligible withrespect to the prompt one, the effect of diffusion has aminimal effect in the total flux. The inclusion of diffusionhas the effect of reducing the ICS flux by a renormaliza-tion factor that changes at most of about 20 −
25% theflux for ICS. Since the numerical calculation of Eq. 6,that includes diffusion, is very time consuming and itsaddition does not change significantly the predictions forthe γ -ray flux, we decide to neglect this process in thecalculation. Therefore, we will use Eq. 14 in our analy-sis. We will discuss in Sec. III B how the inclusion of thediffusion can affect our results. γ rays from bremsstrahlung There is an additional secondary production of γ raysfrom DM that is associated with the Bremsstrahlung pro-cess. This involves e ± , produced from DM particles anni-hilation, interacting with interstellar gas, in the neutral,ionized and molecular forms, and generating photons typ-ically at X-ray and γ -ray energies. The calculation ofthis contribution follows the one in Eq. 14, if the diffu-sion is not taken into account, where the γ -ray power forICS is substituted with the one for Bremsstrahlung. Forthis latter quantity we consider the approximated formin Ref. [59].We show in Fig. 1 the contribution of Bremsstrahlung γ rays to the total DM contribution for the µ + µ − and τ + τ − and b ¯ b channels, assuming for the interstellar gasan average density of 1 cm − . This value is justified bythe density of interstellar gas in the inner few kpc fromthe Galactic center. In particular the distribution of gason the Galactic plane in the inner 3-5 kpc is between 1-3cm − . However, the gas density decreases as an expo-nential function with scale radius of about 0.1-0.2 kpc[57]. For the DM mass and cross section we assume 50GeV and 3 × − cm /s. Bremsstrahlung contributesmostly at energies below 1 GeV to the total flux. Theaddition of this mechanisms does not have a significanteffect since for µ + µ − ( τ + τ − ) it is a factor of about 5 (3)smaller than the ICS one and for b ¯ b it is much smallerthan the prompt emission for most of the energies consid-ered. In addition, the effect that Bremsstrahlung bringsto the total flux is opposite with respect to the addi-tion of diffusion. Therefore, the combination of addingBremsstrahlung emission and the diffusion process in thecalculation has the net effect of producing a difference inthe total flux that is minimal with respect to the casewhere we include prompt and ICS emission only withoutaccounting for diffusion. There is an additional reason toassume that Bremsstrahlung is not contributing signifi-cantly to the GCE. The γ -ray flux for Bremsstrahlungwould be associated with the distribution of the inter-stellar gas distribution that is elongated on the Galactic E [GeV]10 E d N / d E [ G e V / c m / s / s r ] + PromptICS no diffICS with diffBrem no diff ICS+Prompt no diffTOT no diffTOT with diff E [GeV]10 E d N / d E [ G e V / c m / s / s r ] + PromptICS no diffICS with diffBrem no diff ICS+Prompt no diffTOT no diffTOT with diff E [GeV]10 E d N / d E [ G e V / c m / s / s r ] bb PromptICS no diffICS with diffBrem no diff ICS+Prompt no diffTOT no diffTOT with diff
FIG. 1: Flux of γ rays produced for prompt (blue dashed line),Bremsstrahlung (dotted green line) and ICS emission fromDM particles annihilating into µ + µ − , τ + τ − and b ¯ b , from topto bottom panels, for M DM = 50 GeV and (cid:104) σv (cid:105) = 3 × − cm /s. We present two cases for the ICS emission with(orange dot-dashed line) and without (red dotted line) ac-counting for diffusion. We also display the total emission(prompt plus ICS and Bremsstrahlung), with (solid blackline) and without (grey dot-dashed line) diffusion and the casewith ICS, calculated without diffusion, plus prompt emission(dashed brown line). The case with prompt emission andICS, calculated without diffusion, is the model we use in theanalysis. plane. However, there is no evidence that there is anasymmetry on the Galactic plane of the GCE.To summarize, we decide to perform the calculationfor the γ -ray flux from DM by including prompt andICS emissions and without accounting for the diffusionmechanism. The ISRF is modeled as in the model ofRef. [54] for the Galactic center region. B. Antiprotons and positrons flux from darkmatter
1. Dark matter and astrophysical source terms
DM annihilation can induce a primary flux of ¯ p and e + .The source term which denotes the differential produc-tion rate of i = ¯ p, e + per volume, time and energy readsexactly as in Eq. 13 where ( dN i /dE ) f is the spectrum ofantiprotons/ positrons for the annihilation channel. Weagain take ( dN i /dE ) f from Ref. [50]. In addition, thereis an astrophysical antimatter background which orig-inates from the scattering of CR protons and nuclei onthe interstellar matter. The source term for this so-calledsecondary production reads: Q sec i = (cid:88) j,k π (cid:90) dE (cid:48) (cid:18) dσ jk → i dE (cid:19) n k Φ j ( E (cid:48) ) , (15)where Φ j denotes the flux of the progenitor species j ,while n k stands for the number density of the target nu-cleus k in the Galactic disc. We can set j, k = p, Hesince contributions from heavier nuclei are strongly sup-pressed. Furthermore, we will approximate the protonand helium fluxes as spatially constant in the Galacticdisc. This simplification leaves local secondary fluxesvirtually unaffected since the radial dependence of pro-genitor fluxes is effectively absorbed into the propagationparameters.The differential ¯ p production cross sections dσ ij → ¯ p /dE are taken from Ref. [60, 61] and include the full modelingof prompt ¯ p emission as well as displaced ¯ p productionvia hyperon and ¯ p decays (see Refs. [62, 63] for otherrecent cross section parametrizations). Since productioncross sections are only known to a few percent precision,they comprise an important source of systematic error inthe modeling of ¯ p fluxes which needs to be incorporatedin DM searches. These uncertainties and their full corre-lations have been parameterized in Ref. [61] and will beincluded in our analysis as well.In the case of cosmic e + , secondary production con-tributes strongly to the astrophysical background con-tribution below 10 GeV while at higher energies the cu-mulative flux from Galactic pulsar wind nebulae (PWNe)liekely dominates. Since the contribution of PWNe to the e + data is still not well constrained, we are interested inproviding conservative constraints to the DM contribu-tion. Therefore, we use the e + production cross section parameterization of Ref. [64] which yields the lowest sec-ondary flux among the parametrizations in the literature(see Ref. [65]).
2. Antimatter propagation
The propagation of antimatter follows a transportequation analogous to Eq. (9). Besides diffusion and en-ergy losses, we, however, also include reacceleration bymagnetic shock waves as well as annihilation processesin the Galactic disc. Convective winds will be neglectedsince they are not preferred by recent CR analyses (seee.g. [47, 66]).We solve the transport equation within the two-zonediffusion model [67, 68, 69] which assumes that diffusionoccurs homogeneously and isotropically in a cylinder ofradius R and half-height L around the Galactic disc. Thedisc itself is taken to exhibit a thickness 2 h = 0 . n H = 0 . − and n He = 0 . − . With theseassumptions, the transport equation becomes: − K ∆ N i + 2 hδ ( z ) (cid:2) ∂ E ( b disc N i − K EE ∂ E N i ) + Γ ann N i (cid:3) + ∂ E ( b halo N i ) = 2 hδ ( z ) Q sec i + Q prim i . (16)The extension of the disc in vertical direction (z-direction) has been neglected. Processes confined to thedisc were multiplied by 2 hδ ( z ) for proper normalization.The diffusion coefficient K is modeled as a brokenpower law in the rigidity R [70] K = K β η (cid:18) R GV (cid:19) δ (cid:32) (cid:18) RR b (cid:19) ∆ δ/s (cid:33) − s , (17)with power law index δ below the break position R b and δ + ∆ δ above. The parameter s describes the smoothnessof the break. The diffusion break is required to accountfor observed spectral breaks in the proton and nuclearcosmic ray spectra [71, 72] but plays a subleading rolein the energy range accessible to antimatter searches. We also allow for a free scaling of K with the velocity β . While η = 1 in the original two-zone diffusion model,recent studies discovered a significant improvement in thefit to secondary nuclear cosmic rays if η is taken as a freeparameter [66, 74, 75, 76]. Physically, an increase of thediffusion coefficient (negative η ) towards low rigidity ismotivated by wave damping on cosmic rays [77].Reacceleration by Alfv´en waves is modeled as diffusionin momentum space via the term [69] K EE = 43 V a K p δ (4 − δ )(4 − δ ) , (18) A break in the diffusion term can be linked to the transitionfrom diffusion on CR self-generated turbulence at low rigidity todiffusion on external turbulence at high rigidity [73]. where V a stands for the Alfv´en velocity. Energy lossesin the Galactic disc arise from Coulomb interactions,ionization, bremsstrahlung and reacceleration, such that b disc = b coul + b ion + b brems + b reac . We extract b coul , b ion , b brems from [78] and b reac from [69].For e + , we also need to include the energy loss term b halo = b ic + b synch which accounts for inverse Comptonscattering and synchrotron emission in the Galactic haloas described in Sec. II A 2. We use for the ICS calculationthe full Klein Nishina formalism and the ISRF model asin Ref. [54].Annihilation in the Galactic disc is only relevant forantiprotons. The annihilation rate Γ ann is taken fromRef. [79, 80]. Furthermore, we consider inelastic (non-annihilating) scattering of antiprotons with the interstel-lar matter through inclusion of a tertiary source term asin [68].The spatial part of the antiproton transport equationcan be solved analytically. For secondary antiprotons,whose source term is located in the Galactic disc, oneobtains [67, 68] (cid:18) h Γ ann + 2 KL (cid:19) N ¯ p + 2 h∂ E ( b disc N ¯ p − K EE ∂ E N ¯ p )= 2 h ( Q sec¯ p + Q ter¯ p ) , (19)with the tertiary source term as defined in [68]. Thisequation needs to be solved numerically. Since the pri-mary antiproton source term contains an additional spa-tial dependence on the dark matter profile, the solutionfor antiprotons from dark matter requires a Bessel ex-pansion in the radial coordinate. The procedure has beendescribed in full detail in [82, 83].An approximate solution of the transport equation forpositrons was already given in Eq. 10. The latter consid-ers only diffusion and halo energy losses, while neglect-ing reacceleration as well as positron interactions withmatter in the Galactic disc. When determining the localcosmic ray positron flux we include the vertical boundaryof the diffusion zone (while the radial boundary can stillbe neglected). The free Green function given in Eq. 11gets modified and one obtains [84] G ( E e , r ← E s , r s ) = 1 b ( E e )( πλ ) / ∞ (cid:88) n = −∞ ( − n × exp (cid:18) − ( x − x s ) + ( y − y s ) + ( z − z sn ) λ (cid:19) , (20)with r = { x, y, z } and z sn = 2 nL + ( − n z s . (21) The antiproton annihilation cross section was interpolated be-tween the two parameterizations as in [81]. We neglected the radial boundary R which is justified since R (cid:28) L for the propagation configuration we consider in this work. The solution in Eq. 10 with the Green function as definedabove holds for both, primary and secondary positrons.In a next step, one can include reacceleration and discenergy losses for positrons through the pinching methoddescribed in Ref. [85]. We note that reacceleration anddisc losses only affect the low-energy range ( E (cid:46) φ , which is universal amongCR species. In this work we will employ an improvedforce field approximation which additionally allows us toinclude charge breaking effects due to CR drifts. Duringa positive solar polarity phase, positively charged parti-cles access the heliosphere on direct trajectories along thepoles. Negatively charged particles enter by inward driftalong the current sheet which gives rise to additional en-ergy losses [87, 88]. In order to incorporate these effects,the Fisk potential is written as a rigidity-dependendentfunction of the form [89, 90] φ = φ + φ GV R . (22)The second term on the right-hand-side models the in-creased energy loss along the current sheet faced by par-ticles whose charge sign is opposite to the solar polarity.It is taken to vanish for particles with charge sign equalto the polarity (in which case the standard force fieldapproximation is recovered). For the AMS-02 data tak-ing period, which (mostly) refers to a positive polarityphase, we will take φ = 0 for positrons, but non-zero forantiprotons. III. DARK MATTER INTERPRETATION OFTHE GALACTIC CENTER EXCESSA. Dark matter density
One of the main ingredients to calculate γ -ray fluxesfrom DM is its density distribution in the Galaxy thatenters through the geometrical factor ¯ J (see Eq. 14). Weuse the surface brightness data of the GCE reported inRef. [14] and the recent results from the rotation curve ofthe Milky Way from Ref. [49] to estimate the values of theDM density parameters. We employ the results obtainedin this section for the estimation of ¯ J for γ rays but alsofor the calculation of the ¯ p and e + production from DM.We derive the predicted surface brightness from DMcalculating the γ -ray flux for different angular distancefrom the center of the Galaxy, i.e. in Eq. 14 the geomet-rical factor becomes a function of the angular distancefrom the Galactic center. We test the three DM den-sity profiles reported in Sec. II A 1: gNFW, Einasto andBurkert. For the gNFW and Einasto we fit the values of γ and α , respectively. For both profiles we fix r s = 20kpc since the surface brightness data are at small an-gular distances from the Galactic center and r s is thusunconstrained. In fact, we check that by using differentvalues for r s the results do not change. Finally, for theBurkert profile we leave free to vary r s since the slopeis fixed. The values of normalization of the DM den-sity profile cannot be derived with this method since, inthe flux calculation, ρ s is completely degenerate with theannihilation cross section.We find the best-fit values of γ for the gNFW, α for theEinasto and r s for the Burkert profile by fitting the pre-dicted DM flux to the GCE data for the surface bright-ness recently measured in Ref. [14] between 0 ◦ to 20 ◦ from the Galactic center. The result of the fit is thatthe γ parameter for the gNFW profile must be between1 . − . χ ( ˜ χ ) 3 . . γ = 1 . α = 0 .
13 and ˜ χ = 1 .
9. Finally, the Burkertprofile gives a very poor fit with ˜ χ = 12 .
8. The Burk-ert profile gives a flat surface brightness in the inner fewdegrees from the Galactic center where instead the dataare very peaked. In addition, the best-fit for r s is about0.26 kpc that is a too small value if compared to the ob-served DM density in the outer part of the Galaxy (see,e.g., [50]). Therefore, we decide to consider the followingthree cases to bracket the possible uncertainty of the DMdensity profile: gNFW with γ = 1 . . α = 0 .
13. We show the best-fit we obtain with thesethree models in Fig. 2 compared to the GCE data for thesurface brightness obtained with the
Baseline
IEM. Inparticular we can observe that all the three cases providea good fit to the GCE data. DM is able to fit properlythe peaked data in the inner few degrees from the Galac-tic center but also the extended tail beyond 5 ◦ . We alsotest the same analysis using the surface brightness dataobtained in Ref. [14] with other IEMs and we find verysimilar results for γ and α .We derive the normalization ρ s and the scale radius r s of the DM density profile using the results for the localDM density and total DM mass published in Ref. [49].The authors analyze precise circular velocity curve mea-surements of the Milky Way for distances between 5 − ρ (cid:12) = [0 . , .
39] GeV/cm . Instead, the dS d [ M e V / c m / s / s r ] NFW = 1.2NFW = 1.3Einasto, = 0.13GCE data Di Mauro 2020
FIG. 2: Result of the fit to the GCE surface brightness data(black data points) [14] with a DM signal calculated for agNFW profile with γ = 1 . γ =1 . α = 0 . γ / α ) ρ s [GeV/cm ] r s [kpc] ¯ J label ρ (cid:12) = 0 .
300 GeV/cm M = 5 . · M (cid:12) gNFW 1.20 0.416 12.87 111.5 MIN gNFW 1.30 0.314 14.18 155.3Einasto 0.13 0.376 7.25 288.9 ρ (cid:12) = 0 .
345 GeV/cm M = 5 . · M (cid:12) gNFW 1.20 0.587 11.57 166.1gNFW 1.30 0.449 12.67 231.0 MED
Einasto 0.13 0.569 6.35 449.3 ρ (cid:12) = 0 .
390 GeV/cm M = 6 . · M (cid:12) gNFW 1.20 0.851 10.20 246.8gNFW 1.30 0.649 11.20 339.1Einasto 0.13 0.864 5.51 686.7 MAX
TABLE I: This table summarizes the best-fit for the DM den-sity parameters for each case considered in the paper. We listnine cases that result from choosing three different DM den-sity profiles and three measurements for the local DM densityand M from Ref. [49]. We report the value of the slope ( γ for the gNFW and α for Einasto), ρ s , r s and the value of thegeometrical factor ¯ J calculated for an ROI 40 ◦ × ◦ centeredin the Galactic center. DM mass is provided through the quantity M , definedas the mass contained within the radius r such thatthe energy density is 200 times larger than the criticalenergy density of the Universe. M is found to varybetween [5 . , . × M (cid:12) for r = [175 , γ = 1 . ρ (cid:12) = 0 .
30 GeV/cm and M = 5 . × M (cid:12) (see Tab. 2 of [49]) and γ = 1 . ρ (cid:12) = 0 . and M = 6 . × M (cid:12) (see Tab. 3 of[49]). We introduce a scenario that is an average ofthe previous two and defined as ρ (cid:12) = 0 .
345 GeV/cm and M = 5 . × M (cid:12) . We also use the Einastoprofile with α = 0 .
13, that provides a good fit to theGCE data, with the above cited three set of values for ρ (cid:12) and M . We employ the parameters reported forthe previous three set of ρ (cid:12) and M to estimate thevalue of r s and ρ s . In particular these two parametersare found by fixing the local DM density to the values ρ (cid:12) = [0 . , . , . and by integratingthe DM density as (cid:82) r d rρ ( ρ s , r s ) such that M isequal to M = [5 . , . , . × M (cid:12) , respectively foreach ρ (cid:12) value.We report in Tab. I the best-fit values for ρ s and r s thatwe find applying this technique to the three DM densitymodels used in the paper. Since we assume three DMdensity profiles and we consider three possible choices ofthe quantities ρ (cid:12) and M , we end up with nine possi-ble scenarios for the parametrization of ρ . We calculatefor each of the nine cases the value of ¯ J using Eq. 2.As expected with a larger value of the local DM densityalso the value for the geometrical factor is larger. In par-ticular by looking to the gNFW with γ = 1 . J changes from 155 to 339 by varying ρ (cid:12) from 0.30 to 0.39GeV/cm . This increase is proportional to the variationof ρ (cid:12) . Moreover, by changing the DM density profilefrom gNFW with γ = 1 . γ = 1 . J due to the mod-eling of the DM density, and in particular in its localdensity and functional form. These are the cases gNFWwith γ = 1 . ρ (cid:12) = 0 .
300 GeV/cm , labeled as MIN , γ = 1 . ρ (cid:12) = 0 .
345 GeV/cm , named as MED , andEinasto with ρ (cid:12) = 0 .
390 GeV/cm with MAX . The valueof the geometrical factor varies by a factor of 6.2 betweenthe
MIN and the
MAX models.The variation we consider in this paper for ρ encom-passes the systematic on the choice of the DM densityprofile and the local DM density. However, some of theliterature papers find an even larger variation becauseof estimate of the local DM density that is beyond ourrange of 0 . − .
39 GeV/cm . For example, Refs. [92, 93]report values larger than 0 .
40 GeV/cm . Using the re-sults of these latter papers would have the consequenceof providing smaller values of annihilation cross sectionwith respect to the ones reported in this paper. Our re-sults on the systematic of the DM density distributionare similar to the ones obtained recently in Ref. [94] withthe Milky Way rotation curve data. In particular, theirvariation of the DM density parameters produce a sys-tematics on the value of the geometrical factor similar towhat we estimate using the models MIN , MED and
MAX .We only assume annihilating DM because fitting theGCE surface brightness data with this model providesDM density profiles compatible with expectations fromN body simulations. Instead, the calculation of the ge-ometrical factor for decaying DM would be proportional M DM [GeV]10 v [ c m / s ] e + e ++ qqccbbgg FIG. 3: Best-fit for the DM parameters M DM and (cid:104) σv (cid:105) ob-tained by fitting the GCE data in Ref. [14]. The values ofthese data points are reported in Tab. II. The green datapoint labeled with q ¯ q denotes a DM annihilation channel intothe light quarks u, d, s .Channel M DM [GeV] (cid:104) σv (cid:105) [ × − cm /s] χ ( ˜ χ ) e + e − +4 − . +0 . − . .
61 (5 . µ + µ − +11 − . +0 . − . .
12 (5 . τ + τ − . +1 . − . . +0 . − . .
40 (39 . q ¯ q +4 − . +0 . − . .
89 (6 . c ¯ c +3 − . +0 . − . .
11 (7 . b ¯ b +6 − . +0 . − . .
47 (5 . gg +3 − . +0 . − . .
14 (7 . M DM and (cid:104) σv (cid:105) derived by fitting the GCE data obtainedin Ref. [14] with different IEMs. The errors on M DM and (cid:104) σv (cid:105) represent the variation of the best-fit values due to thesystematic on the IEMs. We also display the value of the χ ( ˜ χ ). to ρ . Therefore, in order to fit well the GCE data, valuesaround γ ∼ . γ ∼ B. Fitting the Galactic center excess SED
In this section we fit the GCE SED measured inRef. [14] in order to find the best-fit DM mass and an-nihilation cross section. We use Eq. 14 to calculate the γ -ray flux for the prompt and ICS emission.0 E [GeV]10 E d N d E [ G e V / c m / s / s r ] e + e E [GeV]10 E d N d E [ G e V / c m / s / s r ] bb FIG. 4: Best-fit γ -ray flux obtained with e + e − (top panel) and b ¯ b (bottom panel) annihilation channels (blue dashed line)compared to the GCE data (black data points) reported inRef. [14]. The grey band takes into account the variation inthe GCE data found by performing the analysis with differentanalysis techniques and IEMs.
1. Single channel case
First we assume the simplest scenario with DM parti-cles annihilating into a single channel (i.e. Br = 1). Weconsider the following channels: leptonic ( e + e − , µ + µ − , τ + τ − ), quarks q ¯ q ( q = u, d, s denotes a light quark), c ¯ c , b ¯ b and gluon Gauge bosons gg . All the plots and χ values are found by fitting the GCE data obtained inRef. [14] with the Baseline
IEM. The case with the t ¯ t quark, the Gauge bosons Z Z , W + W − and the Higgsbosons hh provide very poor fits to the GCE flux sincethe masses of these particles are higher than at least 80GeV and the GCE flux peaks at much smaller energies.Therefore, we decide to avoid reporting the results weobtain with these channels.In Tab. II and Fig. 3 we show the results for the bestfit of M DM and (cid:104) σv (cid:105) . The errors represent the variationof the DM parameters derived by fitting the GCE SEDdata obtained with the different IEMs in Ref. [14]. The annihilation channels that provide the best match withthe data, with increasing values of the chi-square ( χ ),are: e + e − , µ + µ − , b ¯ b , q ¯ q , c ¯ c , g ¯ g and τ + τ − . The reducedchi-square ˜ χ = χ /d.o.f. obtained for the quarks chan-nels b ¯ b , c ¯ c , q ¯ q is between 6 and 7 while for the e + e − and µ + µ − ones is about 5.4. The channel τ + τ − , instead, pro-vides a much poorer fit with ˜ χ = 39 .
3. Therefore, thislatter channel alone is not able to explain sufficiently wellthe GCE SED. The cases c ¯ c , q ¯ q and g ¯ g provide very simi-lar results for the DM parameters and goodness of fit. Infact, the intrinsic γ -ray spectrum dN γ /dE is very sim-ilar for these channels (see Fig. 3 of [50]). The resultsobtained for the single channel are similar to the onespublished, for example, in Refs. [9, 10].As shown in Sec. II A the inclusion of diffusion pro-cess and Bremsstrahlung emission in the calculation canslightly affect the results. In particular, considering theDM masses we find from the GCE SED, these two in-gredients would change in opposite directions the pre-dictions, i.e. the inclusion of diffusion (Bremsstrahlung)decreases (increases) the predictions for the γ -ray flux.The variations for the best-fit values of M DM and (cid:104) σv (cid:105) depends on the specific value of the gas density consid-ered in the analysis. Assuming a value of about 1 cm − the changes in the best-fit values for the DM couplingparameters are minimal.We show the results obtained with e + e − and b ¯ b an-nihilation channels in Fig. 4. In particular the flux forthe e + e − channel is dominated by the ICS contributionthat has a peak at about a few GeV. Instead, for the b ¯ b channel the SED is mainly due to the prompt emis-sion. As expected, the peak of the prompt emission forthe b ¯ b channel is at about a factor of 10 smaller energythan the DM mass. The values of ˜ χ are larger than 1for all channels meaning that the fit is not sufficientlygood. For the hadronic channels the reason is that the γ -ray flux has a strong softening above roughly 1/10 ofthe DM mass. Therefore, the γ -ray flux above 10 GeVis much smaller than the GCE data (see bottom panelof Fig. 4 for the b ¯ b channel). Instead, the leptonic chan-nels e + e − and µ + µ − gives a larger contribution above10 GeV thanks to interplay between the ICS and promptemission. However, the SED from DM is systematicallyabove the data between 0 . − . e + e − channel improves significantly with a ∆ χ = 25 and arenormalization of the ISRF of 0.70. The best-fit valuesfor M DM and (cid:104) σv (cid:105) become 28 . . × − cm /s. The improvement with µ + µ − gives ∆ χ = 10, arenormalization of the ICS flux of 1.4, M DM = 56 GeV1and (cid:104) σv (cid:105) = 2 . · − cm /s. Instead, the τ + τ − channelcontinues to provide a poor fit to the GCE SED. The b ¯ b channel fit improves by ∆ χ = 20 with a renormal-ization of the ICS flux of 4.5 but the best-fit values for M DM and (cid:104) σv (cid:105) remain unchanged with respect to Tab. II.Variations of the ISRF density of the order of 30% fromthe one in Ref. [54] are possible considering the currentuncertainties in modeling the ISRF in the center of theMilky Way (see, e.g., the differences between the modein Ref. [54] and [95]).
2. Two and three channels cases
In this section we investigate a more complicated sce-nario where DM particles annihilate into two or three an-nihilation channels. In order to account for these caseswe use a branching ratio Br that multiplies the annihila-tion cross section of the first channel, as in Eq. 14, whilethe second channel is multiplied by 1 − Br . For example,a case with Br = 0 . µ + µ − − b ¯ b case implies that (cid:104) σv (cid:105) is multiplied for 0.7 for the former and 0.3 for thelatter channel as follow: dN γ dE = Br dN µ + µ − dE + (1 − Br ) dN b ¯ b dE (23)The procedure we use to find the DM coupling param-eters is the same applied for the single channel in theprevious section. In Tab. III we show the best-fit valuesfor the DM parameters M DM , (cid:104) σv (cid:105) and Br found by fit-ting the GCE flux data obtained with the Baseline
IEM.Instead in Tab. IV we show the uncertainties for the sameparameters derived when we fit the DM flux to the GCEdata obtained with different IEMs as in Ref. [14]. We donot consider here DM annihilating into q ¯ q and g ¯ g sinceit gives very similar results to c ¯ c (see Tab. II).The DM candidates that provide the largest improve-ment in the goodness of fit with respect to Sec. III B 1are µ + µ − − b ¯ b and τ + τ − − b ¯ b with ∆ χ of 74 and 82, re-spectively. These values of ∆ χ are associated with theadditional parameter Br and they imply 8 . . σ significance for the two channels with respect to the sin-gle one. The DM parameters required to fit the GCEflux data are M DM ∼
50 (35) GeV, (cid:104) σv (cid:105) ∼ × − (1 . × − ) cm /s and Br ∼ . µ + µ − − b ¯ b ( τ + τ − − b ¯ b ) DM candidate. Other cases provide a signifi-cant improvements such as c ¯ c − b ¯ b , e + e − − b ¯ b and e + e − − c ¯ c at the 7 .
7, 5 . σ level. In Fig. 5 we show the best fit weobtain for µ + µ − − b ¯ b , τ + τ − − b ¯ b and c ¯ c − b ¯ b . In partic-ular we see that the two channels provide a better fit tothe GCE flux data because the total contribution of γ -ray from DM cover also the energies between 10-30 GeVwhere the single channel was not able to contribute sig-nificantly. Instead, the channels µ + µ − − τ + τ − , µ + µ − − c ¯ c and τ + τ − − c ¯ c , do not provide any improvement in the fitsince the branching ratio value is 0 or 1, i.e. they providea fit with the same χ of the single channel with µ + µ − or c ¯ c . We also test a possible variation of the ISRF densitythat could change the ICS contribution. We perform a fitto the GCE flux by adding a free parameter for the ICScomponent. We find that the goodness of fit improvessignificantly for the e + e − − c ¯ c , e + e − − b ¯ b and µ + µ − − b ¯ b with a ISRF renormalization with respect to the model inRef.[54] of 0.33, 0.10 and 0.10, respectively. The best-fitvalues found for the ICS renormalization are equivalentof reducing the starlight density in the inner part of theMilky Way. Values of 0 . − . χ we find for these three cases are 3.1, 1.8 and1.8 so the fit improves significantly (∆ χ = 24 , , e + e − − b ¯ b that best fits theGCE SED in Fig. 6. We can see that the fit improves sig-nificantly, with respect to the case with renormalizationequal to 1, because with a fainter ICS flux, the low en-ergy flux is more compatible with the GCE data and theprompt emission for the e + e − channel reproduced verywell the flux above 10 GeV that is difficult to fit in themodels tested before.We finally test whether three annihilation channels im-prove further the fit. We consider all the possible com-binations of the single channels reported before. We donot find any significant improvement with respect to thetwo channel cases. In particular, the DM candidate withthe largest improvements are µ + µ − − τ + τ − − b ¯ b withbest-fit parameters M DM = 40 GeV, (cid:104) σv (cid:105) = 1 . × − cm /s, Br = 0 . Br = 0 . µ + µ − − τ + τ − − b ¯ b for M DM = 40 GeV, (cid:104) σv (cid:105) = 1 . × − cm /s, Br = 0 . Br = 0 .
15. These DM candidates improve the fit by3 . σ and 2 . σ significance with respect to the two chan-nel case. IV. DWARF SPHEROIDAL GALAXIESCONSTRAINTS ON THE GALACTIC CENTEREXCESS
In this section we investigate whether the DM can-didates that explain GCE would generate a detectablesignal in the analysis of data from dSphs. We considerfor this scope the list of 48 dSphs published in [28] andthe best-fit values and errors for the geometrical factorsreported in Tab. 1 and A2. We exclude from the listthe satellites of the Andromeda galaxy. We also test thesample of 41 dSphs used in Ref. [29]. We select all theobjects listed in Tab. 1 of Ref. [29] except for the ones la-beled as “Ambiguous Systems”. We take the best-fit anderrors for the geometrical factors as in Tab. 1 for sources For DM particles annihilating in three channels, Br multipliesthe annihilation cross section for channel 1, Br multiplies theannihilation cross section for channel 2 and 1 − Br − Br mul-tiplies the annihilation cross section for channel 3. Channel 1 Channel 2 M DM (cid:104) σv (cid:105) Br χ ( ˜ χ ) ∆ χ (sign.)[GeV] [10 − cm /s] e + e − µ + µ − . ± .
66 1 . ± .
07 0 . ± .
05 126 . .
37) 18(4 . σ ) e + e − τ + τ − . ± .
58 0 . ± .
01 0 . ± .
03 113 . .
92) 31(5 . σ ) e + e − c ¯ c . ± .
57 0 . ± .
02 0 . ± .
05 112 . .
87) 32(5 . σ ) e + e − b ¯ b . ± .
89 1 . ± .
03 0 . ± .
07 112 . .
89) 32(5 . σ ) µ + µ − τ + τ − . ± .
72 3 . ± .
05 1 . ± .
00 164 . .
66) 0(0 σ ) µ + µ − c ¯ c . ± .
72 3 . ± .
05 1 . ± .
01 164 . .
66) 0(0 σ ) µ + µ − b ¯ b . ± .
92 2 . ± .
14 0 . ± .
05 90 . .
12) 74(8 . σ ) τ + τ − c ¯ c . ± .
27 0 . ± .
01 0 . ± .
04 214 . .
38) 0(0 σ ) τ + τ − b ¯ b . ± .
98 1 . ± .
03 0 . ± .
02 82 . .
83) 82(9 . σ ) c ¯ c b ¯ b . ± .
48 1 . ± .
05 0 . ± .
04 115 . .
97) 61(7 . σ ) e + e − c ¯ c . ± .
55 1 . ± .
18 0 . ± .
03 88 . .
07) 56(7 . σ ) e + e − b ¯ b . ± .
81 2 . ± .
17 0 . ± .
03 51 . .
78) 79(8 . σ ) µ + µ − b ¯ b . ± .
95 3 . ± .
21 0 . ± .
02 50 . .
75) 58(7 . σ )TABLE III: This table reports the best-fit for the DM parameters M DM , (cid:104) σv (cid:105) and Br derived by fitting the GCE data inRef. [14] obtained with the Baseline
IEM. The annihilation cross section multiples Br for channel 1 and (1 − Br ) for channel2 as reported in Eq. 23. We also display the value of the χ ( ˜ χ ) and in last column the difference of χ (significance) betweenthe case of the two channel and the single channel reported in Tab. II. The last three rows represent the results we find if weleave free to vary the ISRF density with a renormalization factor with respect to the model in Ref.[54]. The best-fit values forthis renormalization factor is for the three DM candidates, from top to bottom: 0.33, 0.10 and 0.10.Channel 1 Channel 2 M DM (cid:104) σv (cid:105) Br [GeV] [10 − cm /s] e + e − µ + µ − . +16 . − . . +1 . − . . +0 . − . e + e − τ + τ − . +4 . − . . +0 . − . . +0 . − . e + e − c ¯ c . +6 . − . . +0 . − . . +0 . − . e + e − b ¯ b . +6 . − . . +0 . − . . +0 . − . µ + µ − τ + τ − . +12 . − . . +0 . − . . +0 . − . µ + µ − c ¯ c . +11 . − . . +0 . − . . +0 . − . µ + µ − b ¯ b . +9 . − . . +1 . − . . +0 . − . τ + τ − c ¯ c . +3 . − . . +0 . − . . +0 . − . τ + τ − b ¯ b . +5 . − . . +0 . − . . +0 . − . c ¯ c b ¯ b . +4 . − . . +0 . − . . +0 . − . TABLE IV: Same as Tab. III but for the fit performed on theGCE data obtained in Ref. [14] with different IEMs. There-fore, the errors on the DM parameters are due to the variationin the results obtained by fitting the GCE SED data obtainedin Ref. [14] with a variation of the choice for the interstellaremission. with a measured J factor while for the others we usethe value predicted by the J factor-distance relation inEq. 2 of their paper and assuming an error on log ( J )of 0.6. The differences between the sample of dSphs inthe two above cited references are in the list of objectsand the estimated geometrical factors. Bootes III is notconsidered in Ref. [28] while for Tucana III only upperlimits for the geometrical factor are reported. Thus this latter object is not included in the analysis for the sampleof Ref. [28]. The objects Aquarius II, Carina II, Cetus,Leo T are not listed in Ref. [29] while Segue 2 has thechemical signatures of a dSph, but exhibits a low velocitydispersion and so has not been considered. There are alsodifferences in the best-fit values of the geometrical factorthat, however, is for most of the objects well within the1 σ errors. We find similar results using the two samplesat the 15 −
20% level in the relevant mass range for theDM interpretation of the GCE (see Sec. IV B), i.e. for M DM ∈ [10 , A. Data selection and analysis technique
We select the same exposure time of the GCE analysis[14], i.e. eleven years of Pass 8 data (data processingP8R3). We select SOURCEVETO class events , pass-ing the basic quality filter cuts , and their correspond-ing P8R3 SOURCEVETO V2 response functions, as inRef. [14]. We choose energies between 0.3 to 1000 GeV Mission Elapsed Time (MET): 239557417 s − SOURCEVETO is an event class recently created by the
Fermi -LAT team to maximize the acceptance while minimiz-ing the irreducible cosmic-ray background contamination. Infact, SOURCEVETO class has the same contamination level ofP8R2 ULTRACLEANVETO V6 class while maintaining the ac-ceptance of P8R2 CLEAN V6 class. DATA QUAL > E [GeV]10 E d N d E [ G e V / c m / s / s r ] Br ( + )(1 Br ) ( bb )Total E [GeV]10 E d N d E [ G e V / c m / s / s r ] Br ( + )(1 Br ) ( bb )Total E [GeV]10 E d N d E [ G e V / c m / s / s r ] Br ( cc )(1 Br ) ( bb )Total FIG. 5: Flux of γ rays from DM particle annihilating into twochannels. We show the contribution of both channels and thetotal flux compared to the GCE flux data. E [GeV]10 E d N d E [ G e V / c m / s / s r ] Br ( e + e )(1 Br ) ( bb )Total FIG. 6: Same as Fig. 5 leaving free to vary also a renormal-ization of the ICS contribution. For the case reported in thisfigure a renormalization of the ICS emission of 0.1 is foundfrom the fit. and apply a cut to zenith angles < ◦ between 0.3 to1 GeV and < ◦ above 1 GeV in order to exclude theEarth Limb’s contamination. We model the backgroundwith sources reported in the 10-year Source Catalog(4FGL) which is an extension of the 8-year Source Cata-log (4FGL-DR2) [96] . We also include the latest releasedIEM, namely gll iem v07.fits , and its correspond-ing isotropic template iso P8R3 SOURCEVETO V3 v1.txt .We analyze the 12 ×
12 deg regions of interest (ROI) cen-tered in the dSphs position and choose pixel size of 0 . ×
16 deg in order to include also sources atmost 2 ◦ outside our ROI. We will run the analysis withdifferent choices of some of the assumptions done aboveto see how the results change. In particular, we changethe lower bound of the energy range to 0 . ×
15 deg .The analysis of the DM search in our sample of dSphsfollows the one performed in the past by the Fermi -LATCollaboration on these sources (see, e.g., [30]) or morerecently in the direction of Andromeda and Triangulumgalaxies [1]. We provide a general overview and we re-fer to Refs. [1, 30] for a complete description. We usethe public
Fermipy package (version 0.19.0) to perform abinned analysis with eight bins per energy decade.
Fer- https://arxiv.org/pdf/2005.11208.pdf https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/ A complete discussion about this new IEM can be found at https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/Galactic_Diffuse_Emission_Model_for_the_4FGL_Catalog_Analysis.pdf mipy is a python wrapper of the official Fermitools , forwhich we use version 1.3.8. • ROI optimization . A baseline fit is performed oneach ROI including sources in the 4FGL-DR2 cat-alog, IEM and isotropic template. A refinementof the model is run by relocalizing all the sourcesin the model. We check that the new positionsare compatible with the ones reported in the 4FGLcatalog. Then, we search for new sources with aTest Statistic ( T S ) T S >
25 and distant at least1 ◦ from the center of the ROI. A final fit is thenperformed, where all the SED parameters of thesources, normalization and spectral index of theIEM and normalization of the isotropic componentare free to vary. With this first step we thus havea background model that represents properly the γ -ray emission in the ROI. In fact, in all the ROIsconsidered the residuals found by performing a T S map are at most at the level of √ T S ∼ − • DM SED . The DM source associated with eachdSph is added in the center of the ROI as a pointsource, since their predicted angular extension isfor most of them smaller than the
Fermi -LAT PSF(see, e.g. [28]). A fit is then performed. The SEDfor the dSphs is calculated by performing a fit en-ergy bin by energy bin. Specifically, the SED rungives for each energy bin the value of the likelihoodas a function of the DM energy flux. With theSED information we can thus test every possiblespectrum for the source of interest. • Conversion from energy flux to DM space . SpecificDM candidates are tested. We use the SED infor-mation obtained in step two to calculate, for everyannihilation channel, the likelihood as a functionof annihilation cross section and DM mass values.For a given DM annihilation channel and mass thetheoretical DM SED shape is fixed and for differentvalues of (cid:104) σv (cid:105) we extract the correspondent likeli-hoods from the SED data. • Extracting the
T S for the detection of DM or up-per limits for (cid:104) σv (cid:105) . The DM detection T S is foundby finding the minimum of the likelihood in (cid:104) σv (cid:105) and M DM space and comparing it with the like-lihood of the null hypothesis, i.e. the one of theoptimized ROI fit without the DM emission. Theupper limits of (cid:104) σv (cid:105) are instead calculated in thefollowing way. For a fixed DM mass, we take thelikelihood profile as a function of (cid:104) σv (cid:105) ( L ( (cid:104) σv (cid:105) )).We then can calculate the upper limits for (cid:104) σv (cid:105) The Test Statistic (
T S ) is defined as twice the difference inmaximum log-likelihood between the null hypothesis (i.e., nosource present) and the test hypothesis:
T S = 2(log L test − log L null ) [97]. M DM [GeV]0246810121416 T S + Leo VBootes IITucana IIWillman 1 Reticulum IIStacked analysis95% containment68% containment M DM [GeV]0.02.55.07.510.012.515.017.520.0 T S bb Leo VBootes IITucana IIWillman 1Reticulum II Horologium IILeo IStacked analysis95% containment68% containment
FIG. 7:
T S as a function of mass for the dSphs detected withthe highest significance. We also show the
T S for the jointlikelihood analysis on the dSphs sample and the 68% and95% containment bands for the random direction runs. Weshow the results for the τ + τ − (top panel) and b ¯ b annihilationchannels (bottom panel). by finding the minimum of L ( (cid:104) σv (cid:105) ) and calculatingthe (cid:104) σv (cid:105) that worsens the best-fit likelihood valueby ∆ L = 2 . /
2, which is associated with the one-sided 95% CL upper limits. This is the same pro-cedure used in several other papers where the fre-quentist approach is employed (see, e.g. [29]). Infinding the
T S or the upper limits for (cid:104) σv (cid:105) we addto the Poissonian term of the likelihood a factorthat takes into account the uncertainty on the J factor (see Eq. 3 in [30]) taken from [28]. B. Results for the detection and upper limits for (cid:104) σv (cid:105) In this section we report the results for the search ofDM in the directions of the dSphs in our sample. First,5we calculate the
T S of each individual source. We showin Fig. 7 the objects for which we find the highest detec-tion significance: Leo V, Tucana II, Willman 1, Reticu-lum II, Horologium II and Bootes I. Among the dSphsselected the one detected with the highest
T S is Reticu-lum II with a mass of 300 (40) GeV, (cid:104) σv (cid:105) = 1 . × − (9 × − ) cm /s for the b ¯ b ( τ + τ − ) annihilation channeland detected with a T S ∼
10, which corresponds to ap-value of 2 . × − (4 . × − ) local, i.e. pre-trials,significance of ∼ . σ (2 . σ ) . These T S are below thereference value of 25 that is usually used by the
Fermi -LAT Collaboration to include a source in the catalogs.In order to verify more precisely if our findings are sig-nificant or not, we run the same analysis in 100 randomdirections in each ROI. The analysis pipeline is run ex-actly as before but the dSphs emission is searched in otherdirections where we do not expect to detect any signalfrom DM. These simulations provide thus the expectedsignal in case of the null hypothesis, i.e. that there is noemission from DM in dSphs. In Fig. 7 we show the 68%and 95% containment bands for the
T S for the runs inthe 100 random directions. The
T S profiles found formost of the dSphs are compatible with the results of therandom directions except for Reticulum II, Bootes II andWillman 1. Once we have the likelihood profile for eachdSph as a function of DM mass and annihilation crosssection, we can sum all of them together and get thejoint combined likelihood profile for the entire sample ofdSphs. The result for the
T S as a function of mass forthe joint likelihood analysis is not completely containedinside the 95% containment band of the random directionruns.Since the signal detected from each individual dSphand for the stacked sample does not seem to be signifi-cant, we calculate upper limits for the annihilation crosssection. We display them in Fig. 8 for the b ¯ b and τ + τ − annihilation channels. The 95% CL upper limits are be-low the thermal cross section up to roughly 100 GeV forboth channels. We also display the upper limits obtainedwith the list of dSphs in Ref. [29] and using the geomet-rical factors reported in that publication. The resultsobtained with dSphs in Ref. [29] are similar to the onefound with our reference sample. We also show the 68%and 95% containment bands for the limits obtained in100 random directions. These expected limits in case ofno detection are wider at low mass where the LAT ismore sensitive and could pick up residuals due to faintsources or mismodeling of the IEM. Moreover, the 68%containment band is much narrower than the 95% one,as expected. The limits found for the dSphs are compat-ible with the 95% containment band for both the b ¯ b and τ + τ − annihilation channels. Instead, the observed limits In order to convert the
T S into the p-value and the detectionsignificance, we have considered the analysis in 4800 random di-rections and derived the
T S distribution of the detection of thedSphs. M DM [GeV] v [ c m / s ] + Data, this workAlbert+17 SampleThermal v
95% containment68% containment M DM [GeV] v [ c m / s ] bb Data, this workAlbert+17 SampleThermal v
95% containment68% containment
FIG. 8: 95% CL upper limits for (cid:104) σv (cid:105) for the τ + τ − (top panel)and b ¯ b (bottom pannel) annihilation channels found with thedSphs sample in Ref. [28] (black solid line) and Ref. [29] (Al-bert+17, green dashed line). We also show the 68% (yellowband) and 95% (cyan band) containment band for the limitsobtained in random directions (read the main text for fur-ther details). We report the thermal cross section taken fromRef. [98]. are significantly higher than the 68% containment bandbetween about 50 − b ¯ b and 10 −
200 GeVfor τ + τ − because at these DM masses there is a smallsignal in the joint likelihood analysis as shown in Fig. 7.In Fig. 9 we show the ULs obtained for different as-sumptions of our analysis. In particular we perform theanalysis with the SOURCE IRFs, with a wider ROI of15 ◦ × ◦ , selecting data above 0.5 GeV, and using thedSphs sample from Ref. [29] (Albert2017). The resultsare similar for all the cases reported and in the DM massrange 1-100 GeV that is the relevant one for the DM in-terpretation of the GCE. This implies that our results donot change significantly making different choices of thedata analysis or using a different dSphs sample.Our results for the upper limits with dSphs are similar6 M DM [GeV] v [ c m / s ] bb Pace&Strigari+19Pace&Strigari+19, SourcePace&Strigari+19, 15 × 15 ROIPace&Strigari+19, E > 0.5 GeVAlbert+17Thermal v FIG. 9: 95% CL upper limits for (cid:104) σv (cid:105) for the b ¯ b annihi-lation channel for our baseline analysis (Pace&Strigari+19[28], black solid). We also show the limits obtained with theSOURCE IRFs (dashed blue), with a wider ROI of 15 ◦ × ◦ (red dot-dashed), selecting data above 0.5 GeV (green dot-ted), and using the dSphs sample from Ref. [29] (Albert+17,orange solid). at the 20 −
30% level with recently published in Refs [31,32] where different list of sources and analysis techniqueshave been applied.
C. Combining the Galactic center excess withdSphs limits
If DM is responsible for the GCE, an interesting ques-tion arises about its compatibility with the non detectionof a signal from dSphs. In order to answer this question,we compare the coupling parameters of the DM candi-dates that explain the GCE with the limits found fromdSphs. We test the one/two/three channels cases thatprovide the best fits to the GCE SED: b ¯ b and µ + µ − , τ + τ − − b ¯ b and µ + µ − − τ + τ − − b ¯ b . We take the valuesof the masses, annihilation cross sections and branchingratio from Tabs. II and IV that contain the systematicdue to the choice of the IEM. For the first time in liter-ature the limits for (cid:104) σv (cid:105) for dSphs are calculated assum-ing specific models with two and three annihilation chan-nels. This is done with the same procedure explained inSec. III B but assuming for the intrinsic γ -ray spectrumfrom DM dN γ /dE the specific DM two or three channelbranching ratios.We show the result of this analysis in Fig. 10. The GCEDM candidate obtained with the µ + µ − is below the lim-its, even in the 68% CL level case, which is the strongest.However, we have to stress that in the calculation of the γ rays from the GCE we have included both the ICS andprompt emission while for the flux from dSphs we haveaccounted only for the prompt emission. For DM with a mass of 60 GeV the peak of the emission is at abouta few GeV and it is mainly due to ICS on starlight (seetop panel of Fig. 5). Since the stellar light in dSphs isorders of magnitude smaller than in the Milky Way, theICS contribution is negligible with respect to the promptemission. Instead, the annihilation channels b ¯ b , τ + τ − − b ¯ b and µ + µ − − τ + τ − − b ¯ b are dominated by the prompt γ -ray emission from the b ¯ b annihilation channel. Thus theeffect of the diffusion in the ICS calculation for dSphs,that we do not take into account in our calculation, isnegligible. For all these channels the properties of theDM candidate that explains the GCE in the MED DMmodel is roughly at the 95% CL upper limits of the dSphslimits. This implies a tension at about 2 σ significance.However, considering the variation in (cid:104) σv (cid:105) obtained byconsidering the MIN and MAX models, the GCE inter-pretation of DM is compatible with the 68% CL upperlimits of the dSphs, that implies there is no tension. V. CONSTRAINTS ON DARK MATTER USINGAMS-02 ¯ p DATA
Messengers that have provided tight constraints onDM in the past are ¯ p CRs. It is thus very interesting toinvestigate the compatibility of the DM interpretation ofthe GCE with the newest ¯ p flux data collected in 7 yearsof mission by AMS-02 [35]. This is particularly true sincea tentative DM signal has previously been found in theAMS-02 ¯ p data [46, 99] which was argued to be compat-ible with the GCE [43, 45]. On the other hand, it wasnoted that the significance of the ¯ p excess is drastically re-duced, once uncertainties in the production of secondaryantiprotons [44, 48] and the correlations in the AMS-02systematic errors [47, 100] are properly included.We will perform our ¯ p analysis mostly following theapproach described in [47, 48]. The main aspects shallbriefly be described below. In a first step, the high energybreak in the diffusion coefficient (Eq. (17)) is fixed by a fitto the primary proton, helium, carbon, nitrogen and oxy-gen fluxes fluxes measured by AMS-02 [48]. The breakparameters take the values R b = 275 GV, ∆ δ = 0 . s = 0 . p spectrum in the energy range coveredby data, uncertainties in the break parameters can beneglected for our purposes.The Fisk potential parameter φ for the AMS-02data taking period from a combined fit to the AMS-02[101] and Voyager [102] proton data falls in the range φ = 0 . − .
72 GV. The uncertainty encompasses dif-ferent parameterizations of the interstellar proton flux,while statistical errors are negligible [48]. For the sakeof a conservative approach we adopt the upper value φ = 0 .
72 GV in the following. The diffusion coefficientparameters K , δ , η and the Alfv´en velocity V a are de-termined within a joined fit to the AMS-02 ¯ p and B/Cdata [35].In addition we allow the solar modulation parameter7 M DM [GeV]10 v [ c m / s ] + GCE, Syst. DM densityGCE, Syst. IEMsdSphs ULs, 68% CLdSphs ULs, 95% CLdSphs ULs, 99% CL M DM [GeV]10 v [ c m / s ] bb GCE, Syst. DM densityGCE, Syst. IEMsdSphs ULs, 68% CLdSphs ULs, 95% CLdSphs ULs, 99% CL M DM [GeV]10 v [ c m / s ] B r + + (1 B r )( bb ) GCE, Syst. DM densityGCE, Syst. IEMsdSphs ULs, 68% CLdSphs ULs, 95% CLdSphs ULs, 99% CL M DM [GeV]10 v [ c m / s ] Br + Br + (1 Br Br ) bb GCE, Syst. DM densityGCE, Syst. IEMsdSphs ULs, 68% CLdSphs ULs, 95% CLdSphs ULs, 99% CL
FIG. 10: Comparison between the 95% (red dotted), 90% (blue dot-dashed) and 68% (black dashed) CL upper limits for (cid:104) σv (cid:105) obtained from the analysis of the dSphs in Ref. [28] and the DM candidate that fit the GCE flux data obtained in with differentIEMs (green data point). We also display with a green band the variation in (cid:104) σv (cid:105) due to the modeling of the DM density inthe inner part of the Galaxy (see Tab. I). We display DM annihilating into b ¯ b and µ + µ − , τ + τ − − b ¯ b and µ + µ − − τ + τ − − b ¯ b channels. φ which accounts for charge-breaking effects to float(Eq. (22)). In order to constrain it we also include the¯ p flux ratio between AMS-02 and PAMELA [103] in ourfit since PAMELA was run in a phase of opposite so-lar polarity compared to AMS-02 (see [48] for details).The last remaining propagation parameter, the verticalhalf-height of the diffusive zone L , cannot be determinedwithin our fits due to a well-known degeneracy with thediffusion coefficient (which applies to stable secondaryCRs). Based on an analysis of radioactive CRs it has re-cently been determined as L = 4 . +1 . − . kpc for the propa-gation configuration we are employing (QUAINT modelin [104]). We will, therefore, mostly focus on L = 2 − σ range in the following.However, we will also test values of L down to 1 . σ deviation from the preferred value. Wenote that the value of L significantly affects the DM-induced flux of ¯ p because for larger L more DM annihi- lations occur within the diffusion zone and more spaceis available for propagation to the earth. As a result inorder to have the same ¯ p flux from DM for a larger L , asmaller annihilation cross section is required.The secondary ¯ p production is modeled through thecross section parameterization derived in Ref. [61]. Thelatter was obtained by a comprehensive analysis of ¯ p pro-duction in fixed target and collider physics experiments.A potential asymmetry between ¯ n and ¯ p production aswell as exotic channels including hyperons have beentaken into account. The secondary production is subjectto uncertainties a the level of 5 −
10% which have alsobeen derived in [61]. These consist of a (fully correlated)normalization uncertainty N ¯ p = 1 ± .
06 as well as un-certainties with a finite correlation length due to smoothvariations in the cross section. We will include N ¯ p as a fitparameter and map the remaining uncertainties into the¯ p flux through a covariance matrix (following the proce-8dure described in [44, 61, 100]). Similarly, we will includeuncertainties in the boron production through the covari-ance matrix derived in [61].The AMS-02 ¯ p and B/C data exhibit few-percent-levelprecision over wide rigidity ranges. Unless for low andhigh rigidities systematic errors dominate over statisticalerrors. In this light it is unfortunate that correlationsin the systematic errors have so far not been providedby the AMS collaboration. We will estimate the correla-tions in the AMS-02 systematic errors in the ¯ p and B/Cdata following the approach of Ref. [47]. The dominantsystematics come from uncertainties in the CR absorp-tion cross sections within the detector material which aremodeled within the Glauber-Gribov theory in [47]. Wewill also investigate the sensitivity of our results withrespect to the inclusion of correlations in the AMS data.First, we perform a fit to the AMS-02 ¯ p , B/C data andthe antiproton flux ratio between AMS-02 and PAMELAwithout assuming any DM contribution. The best fit ¯ p and B/C spectra are shown together with the AMS-02data in the top panel of Fig. 11. The goodness of fit is χ = 173 on 143 data points with 6 free parameters ofthe model. Therefore, the result for the reduced χ is1.26 which indicates that the AMS-02 data are consis-tent with pure secondary production within ∼ σ . Givensome residual uncertainty in our modeling of correlationsin the AMS systemtic errors (see above) the secondaryhypothesis is definitely in good shape. We report thebest-fit propagation parameters in Tab. V. The cross sec-tion normalization and solar modulation parameters takevalues N ¯ p = 1 .
09 and φ = 0 .
75 GV at the best fit point.The parameters from our fit take values close to thoseobtained in [47] with previous AMS-02 data sets [33, 105].One striking observation is, however, that the residualsbetween the best-fit model and the newest AMS-02 ¯ p fluxdata in the range R = 10 −
20 GV are practically flat.In this rigidity range, previous analyses [46, 99], basedon a previous AMS-02 data set for ¯ p [33], had identifiedthe ‘antiproton excess’ which had tentatively been inter-preted as a DM signal (potentially compatible with theGCE). While the excess occured at a much smaller sig-nificance ( ∼ σ ) after including the correlations in theAMS-02 systematic errors [47], it remained visible in thedata. We realized that the complete disappearance of theexcess is likely linked to the updated AMS-02 data [35]which are systematically lower by ∼
5% in the rigidityrange R = 10 −
20 GV compared to the previous dataset [33].In the next step, we add a DM contribution with freenormalization (cid:104) σv (cid:105) and mass M DM = 7 − bb and ¯ cc are consid-ered. We note that other two-quark as well as two-gluonfinal states yield a very similar ¯ p spectrum as the ¯ cc -channel. Our fits confirm that the previously found ¯ p ex-cess [46, 99] is completely gone in the new AMS-02 data.There is no longer any preference for a DM contribution K [kpc /Myr] δ η V a [km/s]0.042 0.459 -1.49 52.0TABLE V: Best-fit propagation parameters for L = 4 kpcfrom the combined fit to ¯ p and B/C data (assuming pure sec-ondary production of antiprotons). The best fit propagationparameters for different choices of L are obtained by rescaling K with L/ V a by (cid:112) L/ within the range M DM = 30 −
100 GeV. This statementdoes neither depend on the underlying DM profile noron the size of the diffusion zone L which mostly affectthe normalization of a potential DM signal. The bestfit point including a DM contribution is found in the ¯ bb channel at M DM = 1 . ∼ σ ( ∼ σ ) locally (globally).Hence, we do not find any significant preference for aDM signal in the ¯ p data.We can then use ¯ p to provide constraints on DM an-nihilation. Of particular interest is the DM candidate inthe ¯ bb -channel which is compatible with the GCE SED.Employing the parameters reported in Tab. II, we ob-serve that the latter induces a substantial contributionto the the ¯ p flux. If we keep the propagation parametersfixed, we obtain χ = 238 for the MED
DM density modeland L = 3 kpc compared to χ = 173 without DM. If weallow the propagation, solar modulation and cross sectionnormalization parameters to float, χ is reduced to 217which, however, still amounts to an exclusion by > σ forthe DM contribution. The fit in this case prefers a smaller δ = 0 .
36 and higher V a = 63 km/s in order to compen-sate the DM-induced flux which, however, substantiallydegrades the fit to B/C. The best fit ¯ p and B/C spectraincluding the DM contribution are shown in the lowerpanels of Fig. 11. In the following we wish to investigate,whether the exclusion of the GCE DM canditate is ro-bust with respect to variations of the density profile and L .We, therefore, derived the 95% CL upper limits onthe DM annihilation cross section within the mass range M DM = 7 − L = 1 . − MIN , MED , MAX
DM profiles. For the purpose ofderiving limits we keep the propagation parameters fixedat the values indicated in Tab. V, but fully include theuncertainty in the secondary antiproton production crosssection. We tested for a number of parameter points thatallowing the propagation parameters to float would onlyaffect the 95% CL upper limits at the percent level whichis negligible for our purposes.We start by showing the upper limits we find fixingthe DM density model to
MED and testing different val-ues for L . We report this result in the top panel of Fig. 12again for the b ¯ b -channel and the MED
DM density model.We see that the upper limits increase by a factor ∼ L = 5 kpc and 1.5 kpc. This is because forsmall L a large fraction of the ¯ p created at the Galac-9 E$ [GeV]10 E d N / d E [ G e V / c m / s / s r ] SecondaryAMS-02 2018 R [GV]10 B / C r a t i o SecondaryAMS-02 2018 E$ [GeV]10 E d N / d E [ G e V / c m / s / s r ] SecondaryDMTotalAMS-02 2018 R [GV]10 B / C r a t i o SecondaryAMS-02 2018
FIG. 11: Antiproton flux (left column) and B/C (right column) obtained in our fits compared to the AMS-02 data. The upperpanels show the best fit fluxes for pure secondary production (i.e. no DM contribution). The lower panels show the best fitfluxes if we inject a DM signal compatible with the GCE SED for the b ¯ b annihilation channel (with the parameters reportedin Tab. II). The MED
DM density model and L = 3 kpc are assumed. As can be seen, the fit significantly degrades if the DMsignal is added. tic center escapes through the boundaries of the diffu-sion zone before reaching the earth. The DM candidatethat explains the GCE assuming a b ¯ b annihilation chan-nel is compatible with the ¯ p limits only for L ≤ . L constitutes a 3 σ deviation fromthe value preferred by radioactive CRs derived in [104].We also note that another recent evaluation of boron,beryllium and lithium fluxes within a similar propaga-tion setup found L = 6 . ± L ≤ . e + spectrum.In Fig. 12 we also show how the upper limits changeassuming a different DM density distribution. As ex-pected the MIN
DM density provides weaker limits with respect to
MED and
MAX . However, the limits on (cid:104) σv (cid:105) scalealmost proportionally to the change of ¯ J , i.e. the GCEpreferred cross section changes in the same way as the ¯ p limit. Hence, variations in the DM profile do not recon-cile the DM interpretation of the GCE with ¯ p constraints.We have finally tested, whether our conclusions areaffected by the modeling of correlations in the AMS-02data which we adopted from [47]. For this purposes werecalculated the constraints in the ¯ bb -channel assumingsystematic errors are uncorrelated (a common assump-tion in previous CR analyses). However, we found nosignificant change in the limit around M DM ∼
40 GeVcompared to the case where we include AMS-02 correla-tions.We now turn to DM models with a significant annihi-lation fraction into leptons. These should be subject toweaker ¯ p constraints since the antiproton flux from lep-tonic final states is practically negligible. In fact, the DM0 M DM [GeV]10 v [ c m / s ] GCE, MED, bbp , L = 1.5 kpc p , L = 2 kpc p , L = 3 kpc p , L = 4 kpc p , L = 5 kpc M DM [GeV]10 v [ c m / s ] bb GCE, MAXGCE, MEDGCE, MIN p , L = 4 kpc, MIN p , L = 4 kpc, MED p , L = 4 kpc, MAX FIG. 12: Top Panel: 95% CL upper limits on the DM an-nihilation cross section found by fitting AMS-02 ¯ p data andassuming different sizes for L . In addition we show the best-fitDM parameters we obtain by fitting the GCE. We assume a b ¯ b annihilation channel and the MED
DM density model. Bottompanel: Same as top panel changing the DM density model to
MIN , MED and
MAX . The bands we show for the ¯ p upper limitsinclude the variation in the results changing L from 3 to 5kpc. candidates from Tab. II which annihilate into pure e + e − and µ + µ − are not constrained by ¯ p . These channels willbe constrained by CR e + in the next section.However, ¯ p are sensitive to the two-channel final statesof Tab. IV which are partly leptonic and partly hadronic.In Fig. 13 we show 95% CL upper limits for the e + e − − b ¯ b , e + e − − c ¯ c , µ + µ − − b ¯ b and τ + τ − − b ¯ b channels with thebranching ratios from Tab. IV for the MED
DM profile (e.g. e + e − − b ¯ b refers to 50% annihiliation into e + e − and 50%into b ¯ b ). The limits for the b ¯ b -channel are shown again forcomparison. It can be seen that the GCE preferred anni-hilation cross sections are excluded for all mixed channelswith a hadronic component if L = 3 kpc (lower panel ofFig. 13). Reducing the diffusion halo size to L = 2 kpc M DM [GeV]10 v [ c m / s ] MED, p
95% UL for L = 2 kpc GCE, bb GCE, e + e bb GCE, e + e cc GCE, + bb GCE, + bb p , bbp , e + e bbp , e + e ccp + bbp + bb M DM [GeV]10 v [ c m / s ] MED, p
95% UL for L = 3 kpc FIG. 13: Best-fit values for the DM parameters M DM and (cid:104) σv (cid:105) that we find by fitting the GCE SED. We show the casesthat best-fit the GCE SED from Sec. III B. We also report the95% CL upper limits we obtain from ¯ p flux data for the sameDM candidates. We assume the MED
DM density model and L = 2 ( L = 3) kpc for the plot in the top (bottom) panel. reconciles the GCE candidates in the e + e − − b ¯ b , e + e − − c ¯ c and µ + µ − − b ¯ b channels with the ¯ p constraints (upperpanel of Fig. 13). The constraints on the τ + τ − − b ¯ b channel are somewhat stronger due to the larger branch-ing fraction of 80% into b ¯ b (see Tab. IV). We verified thatthese findings remain valid for different choice of the DMprofile.To summarize, all GCE DM candidates which anni-hilate partly or fully hadronically are in some tensionwith the ¯ p constraints. A small diffusion halo L ≤ L ≤ . b ¯ b -channel appears to be the only possible option to rec-oncile the GCE DM candidates with ¯ p constraints. Aswe noted earlier such a small diffusion halo is compatiblewith the observed ¯ p -flux, but causes strong trouble withcomplementary astrophysical probes, in particular withradio data [108, 109, 110] and observations of radioactive1CRs [104, 106]. VI. CONSTRAINTS ON DARK MATTERUSING e + DATA CR e + measured by AMS-02 have been used in thepast to put severe constraints on the leptonic annihila-tion channels of DM. In Ref. [41], for example, the au-thors have assumed that the astrophysical backgroundwas given by an analytic function that was fitting per-fectly the data. They calculated upper limits for (cid:104) σv (cid:105) adding a DM contribution on top of this backgroundmodel. They used this procedure for the leptonic DMchannels, e ± , µ ± and τ + τ − , for which the e + flux shapeis significantly different from the one of the AMS-02 data.However, the resulting constraints can be too optimistic,i.e. too low, because the astrophysical contribution ismodeled by a function that (by construction) perfectlyfits the data and thus almost no space is left for a DMcontribution. In Ref. [42] the authors have done the morerealistic assumption that the e + flux is given by the fol-lowing astrophysical contributions: the secondary pro-duction of primary CRs interacting with atoms of theinterstellar medium and the cumulative flux of PWNe inthe ATNF catalog. The upper limits that they found arehigher than the ones from Ref. [41] but, for the leptonicchannels, they are below the thermal cross section upto about 60 −
100 GeV. The ATNF catalog has a largeincompleteness for sources farther than a few kpc fromthe Earth [111]. These latter sources would mostly con-tribute to the e + flux data below 100 GeV. This energyrange is relevant for a possible contribution of e + fromDM particles with masses below a few hundreds of GeV.In order to account properly for the flux of e + injectedfrom all Galactic pulsars one should perform simulationsbased on synthetic pulsar models (see, e.g., Ref. [112]).Moreover, the secondary production is affected by sys-tematic due to the modeling of the e ± production crosssections usually taken from the one in Ref. [64]. This lat-ter reference, as well as others on the same topic, tunedthe cross sections for the production of e ± with MonteCarlo event generators or old particle data taken decadesago and affected by large statistical and systematic er-rors.A more realistic estimation of the pulsar contributionto the e + as well as the refinement of the e + productioncross sections relevant is beyond the scope of this paper.Therefore, we decide to make two simplistic assumptionsto derive upper limits on the DM annihilation cross sec-tion with AMS-02 e + data [35]. In the conservative ap-proach we assume that the astrophysical e + backgroundis only given by the secondary production, i.e. there isno PWN contribution. Then, we add the DM flux of e + and we use a χ calculation that penalizes modelsthat overshoot the AMS-02 data points. Specifically, ifthe flux from the secondary production and DM is be-low the AMS-02 data the χ remains unchanged, instead E [GeV]10 E [ G e V / m / s / s r ] LPPLETOT Sec. L = 1.5 kpcSec. L = 4 kpc Sec. L = 6 kpc e + AMS-02
FIG. 14: AMS-02 e + flux data (black data points) fitted inthe optimistic approach with an analytic function (black solidline) given by sum of a LogParabola (LP, grey dotted line)and a power-law with an exponential cutoff (PLE, orangedotted line). We also show the secondary flux of e + calcu-lated using the best-fit propagation parameters in Tab. V and L = 1 . , , . − .
82 GV. if it is above the data it is incremented by the typicalfactor (model − data) / (data error) . We show in Fig. 14the comparison between the secondary production calcu-lated for L = 1 . , , . − .
82 GV and the e + data. We use for this anal-ysis the propagation parameters found in Tab. V and aconservative uncertainty of 0.1 GV on the best-fit valueof the Fisk potential obtained by fitting CR data. TheAMS-02 data below 1 GeV rule out vertical sizes of thediffusive halo smaller than 3 kpc. This provides anotherargument against the small value of L required to recon-cile the hadronic GCE DM candidates with ¯ p constraints(see previous section). We test that the e + constraintson (cid:104) σv (cid:105) are similar for all L > L = 4 kpc in the following.The optimistic approach involves the usage of a smoothanalytic function that is able to fit the AMS-02 data.Then, we add the DM contribution and find as 95% CLupper limit the value of (cid:104) σv (cid:105) that worsens the χ fromthe best fit by 2.71. In calculating the best-fit with DMthe free parameters of the analytic functions are left freeto float. This approach is thus similar to the one usedby Ref. [41]. We use a background model that is givenby the superposition of a LogParabola and a power-lawwith an exponential cutoff. This function fits very wellthe data above 1 GeV, in fact the reduced χ is ˜ χ = 0 . conservative and2 M DM [GeV]10 v [ c m / s ] MED, e +
95% UL for L = 4 kpc GCE, e + e GCE, + GCE, bb GCE, e + e cc GCE, + bb GCE, + bb e + , e + ee + , + e + , bbe + , e + e cce + , + bbe + , + bb M DM [GeV]10 v [ c m / s ] MED, e +
95% UL for L = 4 kpc10 M DM [GeV]10 v [ c m / s ] MED, e +
95% UL for L = 1.5 kpc FIG. 15: Best-fit values for the DM parameters M DM and (cid:104) σv (cid:105) that we find by fitting the GCE SED. We show the casesthat best-fit the GCE SED from Sec. III B. We also reportthe 95% CL upper limits we obtain from e + flux data for thesame DM candidates with the conservative (upper panel) andoptimistic methods (central and bottom panels). We assumethe MED
DM density model and L = 4 kpc for the first twopanels and L = 1 . the optimistic approach are shown in Fig. 15 comparedto the best fit of (cid:104) σv (cid:105) and M DM we obtain by fitting theGCE SED. The constraints we calculate for the chan-nel e + e − − b ¯ b are very similar to the ones obtained for e + e − − c ¯ c . The constraints obtained with the conserva-tive approach are compatible with the GCE best fit forall tested cases. As expected the DM annihilation chan-nel with the strongest (cid:104) σv (cid:105) upper limit is the e + e − one.Instead, the results for the optimistic approach are com-patible with the GCE best fit for most single and mixedchannels except for the ones with full or partial annihila-tion into e + e − . In fact, the GCE candidates annihilatinginto e + e − , e + e − − b ¯ b or e + e − − c ¯ c have a cross sectionone order of magnitude higher than allowed by the opti-mistic e + limits. These conclusions do not change if weemploy a lower value of the vertical size of the diffusionhalo L = 1 . e + given by a refined calculation of the secondary produc-tion, tuned on the newest cross section data, and syn-thetic population of pulsars that account properly forthe PWN flux, the upper limits for (cid:104) σv (cid:105) are expected tobe between the ones obtained with the conservative andthe optimistic approach . Therefore, the tension betweenany GCE DM channel with an e + e − contribution andthe AMS-02 e + data is expected to persist even in such amore complete approach. However, since the optimistic e + constraints even for L = 4 kpc only marginally ruleout the dark matter interpretation of the GCE in the µ + µ − -channel, we expect that a more refined analysiswith proper modeling of uncertainties will reconcile thischannel with e + constraints. VII. CONCLUSIONS
In this paper we have shown that the characteristicsof the GCE make DM particles annihilating into theGalactic halo of the Milky Way a viable interpretationfor explaining the excess. In fact, the GCE spatial mor-phology is energy independent and compatible with aNFW profile with γ ∼ . − .
3. Moreover, the GCE isroughly spherically symmetric and its centroid is locatedvery close to the dynamical center of the Galaxy as ex-pected for DM. The GCE SED around the peak at a fewGeV can be well fitted using a single DM annihilationchannel with light quarks, c ¯ c , b ¯ b or the leptonic channels e + e − , µ + µ − with masses from 20 to 60 GeV and crosssections close to the thermal one. We demonstrated thatthe fit to the GCE SED improves significantly in the en-tire energy range by assuming annihilation into two chan-nels with the best cases that are µ + µ − − b ¯ b , τ + τ − − b ¯ b , e + e − − b ¯ b , e + e − − c ¯ c . We have calculated in the pa-per the relevant coupling parameters (mass, annihilationcross section and branching ratio) for each of these cases.Then, we have searched for a cumulative γ -ray signal in Fermi -LAT data compatible with DM particles annihilat-ing in the direction of dSphs. We have performed a com-3bined likelihood analysis of LAT data above 0.3 GeV inwhich we have fully accounted for the uncertainty on theDM density using the information published in Ref. [28]for 48 dSphs. Since we did not find any significant signalwe put upper limits for (cid:104) σv (cid:105) that are below the thermalcross section up to almost 100 GeV for the b ¯ b annihilationchannel. For the first time we tested in a dSphs analysisDM candidates annihilating into two and three channelsfollowing the best-fit cases from the fit to the GCE SED.The upper limits on (cid:104) σv (cid:105) are compatible with the DMinterpretation of the GCE considering the uncertaintiespresent in the DM density distribution.Following a multimessenger approach we have searchedfor a possible DM signal also using the recently released7 years ¯ p and e + AMS-02 flux data. These are amongthe rarest CRs in the Galaxy and have been widely usedin the past as promising cosmic particles for the indirectsearch for DM. First, we analyzed ¯ p data accounting forthe uncertainties in the CR propagation, uncertaintiesin the ¯ p production cross section and the correlation be-tween AMS-02 data points. Since we did not find anysignificant preference for a DM contribution we put up-per limits for (cid:104) σv (cid:105) . The ¯ p constraints exclude all GCEDM candidates reported above with hadronic or semi-hadronic final states unless if the vertical size of the dif-fusive halo is L < − σ below the best fit value obtainedin Ref. [104] using the latest AMS-02 data on radioactiveCRs (see also [106]). Moreover, these small values for L are in tension with complementary astrophysical probes,in particular with radio data [108, 109, 110]. We alsoshowed that variations of the DM density profile cannotreconcile the GCE DM interpretation with ¯ p constraints.Instead, pure leptonic channels are compatible with the¯ p upper limits regardless of the value of L and the as-sumed DM density. Finally, we have calculated upperlimits for a DM contribution from the e + spectrum fol-lowing a conservative approach where only secondary e + were included as background and an optimistic one wherethe e + background is modeled by an analytic function (inorder to also include a potential pulsar contribution). Incase of the conservative approach e + do not provide anyfurther constraints on the DM interpretation of the GCE.In the optimistic approach all mentioned GCE DM candi-dates which annihilate purely or partially into e + e − areruled out.To conclude DM particles annihilating into µ + µ − witha mass of about 60 GeV and a cross section of 4 × − cm /s, which is close to the thermal one, could fit theGCE spectrum. At the same time they are compatiblewith observations of dwarf spheroidal galaxies and would produce a flux of ¯ p and e + compatible with the upperlimits calculated with the latest AMS-02 data. All otherDM annihilation channels we investigated for the GCEare in tension with CR data once we include the latestconstraints on the size of the CR diffusion zone. In par-ticular, the two-channel final state µ + µ − − b ¯ b ( τ + τ − − b ¯ b )with M DM ∼
50 (35), (cid:104) σv (cid:105) ∼ × − ( ∼ . × − )cm /s and Br ∼ . µ + µ − channel, butis compatible with the ¯ p upper limits only for an unfa-vorably small diffusion zone. Acknowledgments
MDM research is supported by Fellini - Fellow-ship for Innovation at INFN, funded by the EuropeanUnion’s Horizon 2020 research programme under theMarie Sk(cid:32)lodowska-Curie Cofund Action, grant agreementno. 754496. MWW acknowledges support by the SwedishResearch Council (Contract No. 638-2013-8993) and bythe Department of Physics of the University of Texas atAustin. MDM acknowledges support by the NASA FermiGuest Investigator Program Cycle 12 through the FermiProgram N. 121119 (P.I. MDM). The authors thankRegina Caputo, Judith L. Racusin, Miguel A. Sanchez-Conde, Michael Gustafsson for providing us commentson the paper.The
Fermi
LAT Collaboration acknowledges generousongoing support from a number of agencies and insti-tutes that have supported both the development and theoperation of the LAT as well as scientific data analysis.These include the National Aeronautics and Space Ad-ministration and the Department of Energy in the UnitedStates, the Commissariat´a l’Energie Atomique and theCentre National de la Recherche Scientifique / InstitutNational de Physique Nucl´eaire et de Physique des Par-ticules in France, the Agenzia Spaziale Italiana and theIstituto Nazionale di Fisica Nucleare in Italy, the Min-istry of Education, Culture, Sports, Science and Technol-ogy (MEXT), High Energy Accelerator Research Organi-zation (KEK) and Japan Aerospace Exploration Agency(JAXA) in Japan, and the K. A. Wallenberg Founda-tion, the Swedish Research Council and the Swedish Na-tional Space Board in Sweden. Additional support forscience analysis during the operations phase is gratefullyacknowledged from the Istituto Nazionale di Astrofisicain Italy and the Centre National d’Etudes Spatiales inFrance. This work performed in part under DOE Con-tract DE- AC02-76SF00515. [1] M. Di Mauro, X. Hou, C. Eckner, G. Zaharijas, andE. Charles, Phys. Rev. D , 123027 (2019), 1904.10977.[2] L. Goodenough and D. Hooper (2009), 0910.2998. [3] D. Hooper and L. Goodenough, Phys. Lett. B697 , 412(2011), 1010.2752.[4] A. Boyarsky, D. Malyshev, and O. Ruchayskiy, Phys. Lett.
B705 , 165 (2011), 1012.5839.[5] D. Hooper and T. Linden, Phys. Rev.
D84 , 123005(2011), 1110.0006.[6] K. N. Abazajian and M. Kaplinghat, Phys.Rev.
D86 , 083511 (2012), [Erratum: Phys.Rev.D87,129902(2013)], 1207.6047.[7] C. Gordon and O. Macias, Phys. Rev.
D88 , 083521(2013), [Erratum: Phys. Rev.D89,no.4,049901(2014)],1306.5725.[8] K. N. Abazajian, N. Canac, S. Horiuchi, andM. Kaplinghat, Phys. Rev.
D90 , 023526 (2014),1402.4090.[9] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden,S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, Phys.Dark Univ. , 1 (2016), 1402.6703.[10] F. Calore, I. Cholis, C. McCabe, and C. Weniger, Phys.Rev. D91 , 063003 (2015), 1411.4647.[11] F. Calore, I. Cholis, and C. Weniger, JCAP , 038(2015), 1409.0042.[12] M. Ajello et al. (Fermi-LAT), Astrophys. J. , 44(2016), 1511.02938.[13] M. Ackermann et al. (Fermi-LAT), Astrophys. J. ,43 (2017), 1704.03910.[14] M. Di Mauro, Submitted to Phys. Rev. D (2021), ArXiv:2101.04694.[15] R. Bartels, S. Krishnamurthy, and C. Weniger, Phys.Rev. Lett. , 051102 (2016), 1506.05104.[16] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer,and W. Xue, Phys. Rev. Lett. , 051103 (2016),1506.05124.[17] O. Macias, C. Gordon, R. M. Crocker, B. Coleman,D. Paterson, S. Horiuchi, and M. Pohl, Nat. Astron. , 387 (2018), 1611.06644.[18] R. Bartels, E. Storm, C. Weniger, and F. Calore, Nat.Astron. , 819 (2018), 1711.04778.[19] R. K. Leane and T. R. Slatyer, Phys. Rev. Lett. ,241101 (2019), 1904.08430.[20] L. J. Chang, S. Mishra-Sharma, M. Lisanti,M. Buschmann, N. L. Rodd, and B. R. Safdi (2019),1908.10874.[21] Y.-M. Zhong, S. D. McDermott, I. Cholis, and P. J. Fox(2019), 1911.12369.[22] S. Abdollahi et al. (Fermi-LAT), Astrophys. J. Suppl. , 33 (2020), 1902.10045.[23] E. Carlson and S. Profumo, Phy. Rev. D , 023015(2014), 1405.7685.[24] J. Petrovi´c, P. D. Serpico, and G. Zaharijaˇs, JCAP , 052 (2014), 1405.7928.[25] D. Gaggero, M. Taoso, A. Urbano, M. Valli, and P. Ul-lio, JCAP , 056 (2015), 1507.06129.[26] N. Aghanim et al. (Planck) (2018), 1807.06209.[27] A. A. Abdo, M. Ackermann, M. Ajello, W. B. Atwood,et al., The Astrophysical Journal , 147 (2010), URL https://doi.org/10.1088/0004-637x/712/1/147 .[28] A. B. Pace and L. E. Strigari, MNRAS , 3480(2019), 1802.06811.[29] A. Albert et al. (Fermi-LAT, DES), Astrophys. J. ,110 (2017), 1611.03184.[30] M. Ackermann et al. (Fermi-LAT), Phys. Rev. Lett. , 231301 (2015), 1503.02641.[31] F. Calore, P. D. Serpico, and B. Zaldivar, JCAP , 029(2018), 1803.05508.[32] S. Hoof, A. Geringer-Sameth, and R. Trotta, JCAP ,012 (2020), 1812.06986. [33] M. Aguilar, L. Ali Cavasonza, B. Alpat, G. Ambrosi,et al. (AMS Collaboration), Phys. Rev. Lett. ,091103 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.091103 .[34] M. Aguilar, L. Ali Cavasonza, G. Ambrosi, L. Ar-ruda, et al. (AMS Collaboration), Phys. Rev. Lett. ,041102 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.041102 .[35] M. Aguilar, L. Ali Cavasonza, G. Ambrosi, L. Arruda,N. Attig, et al., Physics Reports (2020), ISSN 0370-1573, URL .[36] M. Di Mauro, F. Donato, N. Fornengo, et al., JCAP , 006 (2014), 1402.0321.[37] M. Di Mauro, S. Manconi, and F. Donato, Phys. Rev.D , 123015 (2019), 1903.05647.[38] S. Manconi, M. Di Mauro, and F. Donato, Phys. Rev.D , 023015 (2020), 2001.09985.[39] A. U. Abeysekara et al. (HAWC), Science , 911(2017), 1711.06223.[40] M. Di Mauro, S. Manconi, and F. Donato, Phys. Rev.D , 103035 (2020), 1908.03216.[41] L. Bergstrom, T. Bringmann, I. Cholis, D. Hooper,and C. Weniger, Phys. Rev. Lett. , 171101 (2013),1306.3983.[42] M. Di Mauro, F. Donato, N. Fornengo, and A. Vittino,JCAP , 031 (2016), 1507.07001.[43] A. Cuoco, J. Heisig, M. Korsmeier, and M. Kr¨amer,JCAP , 053 (2017), 1704.08258.[44] A. Cuoco, J. Heisig, L. Klamt, M. Korsmeier,and M. Kr¨amer, Phys. Rev. D , 103014 (2019),1903.01472.[45] I. Cholis, T. Linden, and D. Hooper, Phys. Rev. D ,103026 (2019), 1903.02549.[46] M.-Y. Cui, Q. Yuan, Y.-L. S. Tsai, and Y.-Z. Fan, Phys.Rev. Lett. , 191101 (2017), 1610.03840.[47] J. Heisig, M. Korsmeier, and M. W. Winkler, Phys. Rev.Res. , 043017 (2020), 2005.04237.[48] A. Reinert and M. W. Winkler, JCAP , 055 (2018),1712.00002.[49] P. F. de Salas, K. Malhan, K. Freese, K. Hattori, andM. Valluri, JCAP , 037 (2019), 1906.06133.[50] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi,M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Stru-mia, JCAP , 051 (2011), [Erratum: JCAP 10, E01(2012)], 1012.4515.[51] J. F. Navarro, C. S. Frenk, and S. D. M. White, ApJ , 493 (1997), astro-ph/9611107.[52] J. Einasto, Trudy Astrofizicheskogo Instituta Alma-Ata , 87 (1965).[53] A. Burkert, ApJL , L25 (1995), astro-ph/9504041.[54] T. A. Porter, I. V. Moskalenko, A. W. Strong, E. Or-lando, and L. Bouchet, The Astrophysical Journal ,400 (2008), URL https://doi.org/10.1086%2F589615 .[55] S. Blanchet and J. Lavalle, JCAP , 021 (2012),1207.2476.[56] M. Di Mauro and F. Donato, Phys. Rev. D , 123001(2015), 1501.05316.[57] C. Evoli, D. Gaggero, A. Vittino, G. Di Bernardo,M. Di Mauro, A. Ligorini, P. Ullio, and D. Grasso,JCAP , 015 (2017), 1607.07886.[58] M. Cirelli and P. Panci, Nucl. Phys. B , 399 (2009),0904.3830.[59] M. G. Baring, D. C. Ellison, S. P. Reynolds, I. A. Grenier, and P. Goret, ApJ , 311 (1999), astro-ph/9810158.[60] R. Kappl and M. W. Winkler, JCAP , 051 (2014),1408.0299.[61] M. W. Winkler, JCAP , 048 (2017), 1701.04866.[62] M. di Mauro, F. Donato, A. Goudelis, and P. D. Serpico,Phys. Rev. D , 085017 (2014), [Erratum: Phys.Rev.D98, 049901 (2018)], 1408.0288.[63] M. Korsmeier, F. Donato, and M. Di Mauro, Phys. Rev.D , 103019 (2018), 1802.03030.[64] T. Kamae, N. Karlsson, T. Mizuno, T. Abe, and T. Koi,Astrophys. J. , 692 (2006), [Erratum: Astrophys.J.662, 779 (2007)], astro-ph/0605581.[65] T. Delahaye, F. Donato, N. Fornengo, J. Lavalle,R. Lineros, P. Salati, and R. Taillet, Astron. Astrophys. , 821 (2009), 0809.5268.[66] N. Weinrich, Y. G´enolini, M. Boudaud, L. Derome,and D. Maurin, Astron. Astrophys. , A131 (2020),2002.11406.[67] D. Maurin, F. Donato, R. Taillet, and P. Salati, Astro-phys. J. , 585 (2001), astro-ph/0101231.[68] F. Donato, D. Maurin, P. Salati, A. Barrau, G. Boudoul,and R. Taillet, Astrophys. J. , 172 (2001), astro-ph/0103150.[69] D. Maurin, R. Taillet, F. Donato, P. Salati, A. Barrau,and G. Boudoul (2002), astro-ph/0212111.[70] Y. G´enolini et al., Phys. Rev. Lett. , 241101 (2017),1706.09812.[71] M. Aguilar et al. (AMS), Phys. Rev. Lett. , 171103(2015).[72] M. Aguilar et al. (AMS), Phys. Rev. Lett. , 211101(2015).[73] P. Blasi, E. Amato, and P. D. Serpico, Phys. Rev. Lett. , 061101 (2012), 1207.3706.[74] G. Di Bernardo, C. Evoli, D. Gaggero, D. Grasso,and L. Maccione, Astropart. Phys. , 274 (2010),0909.4548.[75] D. Maurin, A. Putze, and L. Derome, Astron. Astro-phys. , A67 (2010), 1001.0553.[76] Y. G´enolini et al., Phys. Rev. D , 123028 (2019),1904.08917.[77] V. Ptuskin, I. V. Moskalenko, F. Jones, A. Strong, andV. Zirakashvili, Astrophys. J. , 902 (2006), astro-ph/0510335.[78] A. Strong and I. Moskalenko, Astrophys. J. , 212(1998), astro-ph/9807150.[79] R. Protheroe, Astrophys. J. , 387 (1981).[80] L. Tan and L. Ng, J. Phys. G , 227 (1983).[81] R. Kappl and M. W. Winkler, Phys. Rev. D , 123522(2012), 1110.4376.[82] A. Barrau, G. Boudoul, F. Donato, D. Maurin, P. Salati,and R. Taillet, Astron. Astrophys. , 676 (2002),astro-ph/0112486.[83] F. Donato, N. Fornengo, D. Maurin, and P. Salati, Phys.Rev. D , 063501 (2004), astro-ph/0306207.[84] E. A. Baltz and J. Edsjo, Phys. Rev. D , 023511(1998), astro-ph/9808243.[85] M. Boudaud, E. Bueno, S. Caroff, Y. Genolini,V. Poulin, V. Poireau, A. Putze, S. Rosier, P. Salati,and M. Vecchi, Astron. Astrophys. , A17 (2017),1612.03924.[86] L. Gleeson and W. Axford, Astrophys. J. , 1011(1968).[87] J. Kota, International Cosmic Ray Conference , 13 (1979).[88] J. R. Jokipii and B. Thomas, Astrophys. J. , 1115(1981).[89] I. Cholis, D. Hooper, and T. Linden, Phys. Rev. D ,043016 (2016), 1511.01507.[90] I. Cholis, D. Hooper, and T. Linden (2020), 2007.00669.[91] A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness,ApJ , 120 (2019), 1810.09466.[92] E. V. Karukes, M. Benito, F. Iocco, R. Trotta,and A. Geringer-Sameth, JCAP , 046 (2019),1901.02463.[93] M. Benito, F. Iocco, and A. Cuoco (2020), 2009.13523.[94] M. Benito, A. Cuoco, and F. Iocco, JCAP , 033(2019), 1901.02460.[95] S. Vernetto and P. Lipari, Phys. Rev. D , 063009(2016), 1608.01587.[96] S. Abdollahi, F. Acero, M. Ackermann, et al., TheAstrophysical Journal Supplement Series , 33(2020), URL https://doi.org/10.3847%2F1538-4365%2Fab6bcb .[97] J. R. Mattox, D. L. Bertsch, J. Chiang, B. L. Dingus,S. W. Digel, J. A. Esposito, J. M. Fierro, R. C. Hart-man, S. D. Hunter, G. Kanbach, et al., ApJ , 396(1996).[98] G. Steigman, B. Dasgupta, and J. F. Beacom, Phys.Rev. D , 023506 (2012), 1204.3622.[99] A. Cuoco, M. Kr¨amer, and M. Korsmeier, Phys. Rev.Lett. , 191102 (2017), 1610.03071.[100] M. Boudaud, Y. G´enolini, L. Derome, J. Lavalle,D. Maurin, P. Salati, and P. D. Serpico, Phys. Rev.Res. , 023022 (2020), 1906.07119.[101] M. Aguilar, D. Aisa, B. Alpat, et al. (AMS Collab-oration), Phys. Rev. Lett. , 171103 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.171103 .[102] E. C. Stone, A. C. Cummings, F. B. McDon-ald, B. C. Heikkila, N. Lal, and W. R. Web-ber, Science , 150 (2013), ISSN 0036-8075,https://science.sciencemag.org/content/341/6142/150.full.pdf,URL https://science.sciencemag.org/content/341/6142/150 .[103] O. Adriani et al., JETP Lett. , 621 (2013).[104] N. Weinrich, M. Boudaud, L. Derome, Y. Genolini,J. Lavalle, D. Maurin, P. Salati, P. Serpico, andG. Weymann-Despres, Astron. Astrophys. , A74(2020), 2004.00441.[105] M. Aguilar, L. Ali Cavasonza, G. Ambrosi, L. Ar-ruda, et al. (AMS Collaboration), Phys. Rev. Lett. ,231102 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.231102 .[106] P. d. l. T. Luque, M. N. Mazziotta, F. Loparco,F. Gargano, and D. Serini (2021), 2101.01547.[107] M. Ackermann et al., ApJ , 3 (2012), 1202.4039.[108] T. Bringmann, F. Donato, and R. A. Lineros, JCAP ,049 (2012), 1106.4821.[109] G. Di Bernardo, C. Evoli, D. Gaggero, D. Grasso, andL. Maccione, JCAP , 036 (2013), 1210.4546.[110] E. Orlando and A. Strong, Mon. Not. Roy. Astron. Soc. , 2127 (2013), 1309.2947.[111] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs,Astron. J. , 1993 (2005), astro-ph/0412641.[112] C.-A. Faucher-Giguere and V. M. Kaspi, Astrophys. J.643