A Robust Excess in the Cosmic-Ray Antiproton Spectrum: Implications for Annihilating Dark Matter

TL;DR

This study confirms a significant excess of cosmic-ray antiprotons around 10-20 GeV, which could be explained by annihilating dark matter particles of about 64-88 GeV, with implications for understanding dark matter and gamma-ray excesses.

Contribution

The paper thoroughly revisits uncertainties in cosmic-ray propagation and solar modulation, confirming the antiproton excess and its potential link to dark matter annihilation.

Findings

Confirmed a 4.7σ antiproton excess after accounting for uncertainties.

Identified dark matter particle mass range of 64-88 GeV consistent with the excess.

The dark matter models also explain the Galactic Center gamma-ray excess.

Abstract

An excess of $\sim$10-20 GeV cosmic-ray antiprotons has been identified in the spectrum reported by the AMS-02 Collaboration. The systematic uncertainties associated with this signal, however, have made it difficult to interpret these results. In this paper, we revisit the uncertainties associated with the time, charge and energy-dependent effects of solar modulation, the antiproton production cross section, and interstellar cosmic-ray propagation. After accounting for these uncertainties, we confirm the presence of a 4.7$\sigma$ antiproton excess, consistent with that arising from a $m_{\chi} \approx 64-88$ GeV dark matter particle annihilating to $b\bar{b}$ with a cross section of $\sigma v \simeq (0.8-5.2) \times 10^{-26}$ cm$^{3}$/s. If we allow for the stochastic acceleration of secondary antiprotons in supernova remnants, the data continues to favor a similar range of dark matter models ($m_{\chi}\approx 46-94$ GeV, $\sigma v \approx (0.7-3.8)\times 10^{-26}$ cm$^3/$s) with a significance of 3.3$\sigma$. The same range of dark matter models that are favored to explain the antiproton excess can also accommodate the excess of GeV-scale gamma rays observed from the Galactic Center.

Tables

Table 1. Table 1: The parameters for the cosmic-ray injection and propagation models used in this study. Each of these models provides a good overall fit to the observed cosmic-ray proton, helium, carbon, and boron-to-carbon ratio spectra up to 200 GV  Cholis  et al. ( 2016 , 2017 ) . These ISM parameters are not picked in advance to fit the p ¯ / p ¯ 𝑝 𝑝 \bar{p}/p ratio.

ISM Model	$\delta$	$z_L\left(\mathrm{kpc}\right)$	$D_0\times {10}^{28}$ (cm2/s)	$v_A$ (km/s)	$d{v}_c/d\left|z\right|$ (km/s/kpc)	${\alpha }_1$	${\alpha }_2$	$R_{\mathrm{br}}$ (GV)
I	0.40	5.6	4.85	24.0	1.0	1.88	2.38	11.7
II	0.50	6.0	3.10	23.0	9.0	1.88	2.45	11.7
III	0.40	3.0	2.67	22.0	3.0	1.87	2.41	11.7

Table 2. Table 2: The values of the dark matter mass and annihilation cross section favored by the AMS-02 antiproton-to-proton ratio, for the case of annihilations to b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} and for each of the cosmic-ray injection and transport models listed in Table  1 . In the top three rows, we have not included any contribution from the acceleration of secondary antiprotons. In the middle three rows, secondary acceleration is included with values of K B subscript 𝐾 𝐵 K_{B} and n gas subscript 𝑛 gas n_{\rm gas} chosen to reproduce the observed boron-to-carbon ratio. In the bottom three rows, secondary acceleration is included with a freely floating value of n gas subscript 𝑛 gas n_{\rm gas} . The rightmost column indicates the statistical preference for a contribution from annihilating dark matter in each case. We remind the reader that we have marginalized over the parameters associated with the antiproton production cross section and solar modulation at each point in the fit (see Sec.  II ) and thus conclude that the excess is statistically significant, even in light of these systematic uncertainties.

ISM Model	$K_B$	$n_{\mathrm{gas}}\left({\mathrm{cm}}^{-3}\right)$	$m_{\chi }\left(\mathrm{GeV}\right)$		$\sigma v\left({\mathrm{cm}}^3/s\right)$	Statistical Preference
I	–	–	$78.3\pm 4.9$		$\left(1.18\pm 0.18\right)\times {10}^{-26}$	4.7$\sigma$
II	–	–	$71.0\pm 3.3$		$\left(2.37\pm 0.28\right)\times {10}^{-26}$	7.7$\sigma$
III	–	–	$74.7\pm 3.8$		$\left(4.17\pm 0.53\right)\times {10}^{-26}$	7.4$\sigma$
I	3.05	2.0 (Fixed)	$81.9\pm 6.1$		$\left(1.08\pm 0.19\right)\times {10}^{-26}$	3.3$\sigma$
II	5.2	2.0 (Fixed)	$57.0\pm 3.1$		$\left(1.63\pm 0.08\right)\times {10}^{-26}$	5.1$\sigma$
III	3.7	2.0 (Fixed)	$51.9\pm 2.8$		$\left(3.05\pm 0.14\right)\times {10}^{-26}$	5.6$\sigma$
I	6.1	0.39 (Float)	$78.1\pm 5.5$		$\left(1.30\pm 0.17\right)\times {10}^{-26}$	3.4$\sigma$
II	10.4	1.28 (Float)	$61.4\pm 2.5$		$\left(1.67\pm 0.10\right)\times {10}^{-26}$	4.6$\sigma$
III	7.4	1.57 (Float)	$52.3\pm 3.3$		$\left(3.20\pm 0.30\right)\times {10}^{-26}$	5.1$\sigma$

Table 3. Table 3: As in Table  2 , but for dark matter that annihilates to b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} , light quarks ( u ​ u ¯ 𝑢 ¯ 𝑢 u\bar{u} , d ​ d ¯ 𝑑 ¯ 𝑑 d\bar{d} ) or W + ​ W − superscript 𝑊 superscript 𝑊 W^{+}W^{-} , and for the case of ISM Model I.

Channel	$m_{\chi }$(GeV)	$\sigma v\left({\mathrm{cm}}^3/s\right)$	Statistical Preference
$b\overline{b}$	$78.3\pm 4.9$	$\left(1.18\pm 0.18\right)\times {10}^{-26}$	4.7$\sigma$
$u\overline{u}$, $d\overline{d}$	$38.0\pm 6.6$	$\left(3.12\pm 0.96\right)\times {10}^{-27}$	3.3$\sigma$
$W^{+}{W}^{-}$	$84.1\pm 2.2$	$\left(3.02\pm 0.42\right)\times {10}^{-26}$	3.6$\sigma$

Equations

$$dN/dR\propto \{\begin{array}{l}{R}^{-{\alpha }_1},\text{for}R<R_{\mathrm{br}}\\ R^{-{\alpha }_2},\text{for}R>R_{\mathrm{br}},\end{array}$$

$$D_{xx}(R)=\beta D_0(R/4\mathrm{GV}{)}^{\delta },$$

$$D_{pp}\propto \frac{R^2{v}_A^2}{D_{xx}(R)},$$

$$v_c=\frac{d{v}_c}{d|z|}|z|.$$

$$\displaystyle N_{CS}(E_{\textrm{kin}})=a+b\,\ln\bigg{(}\frac{E_{\textrm{kin}}}{\rm GeV}\bigg{)}+c\,\bigg{[}\ln\bigg{(}\frac{E_{\textrm{kin}}}{\rm GeV}\bigg{)}\bigg{]}^{2}.$$

$$\frac{d{N}^{\oplus }}{d{E}_{\text{kin}}}(E_{\text{kin}})$$

$$\Phi (R,t,q)$$

$$\rho (r)=\frac{{\rho }_0}{(r/r_s)(1+r/r_s{)}^2}.$$

$${\Gamma }_i(E_{\mathrm{kin}})={\sigma }_i^{\mathrm{inelastic}}\beta c{n}_{\mathrm{gas}}+\frac{1}{E_{\mathrm{kin}}{\tau }_i^{dec}},$$

$$D_i^{\pm }(E)$$

Full text

I Introduction

Measurements of antimatter in the cosmic-ray spectrum have long been used to advance our understanding of high-energy phenomena in the Galaxy Bergstrom et al. (1999); Hooper et al. (2004); Profumo and Ullio (2004); Bringmann and Salati (2007); Pato et al. (2010) . To this end, AMS-02 has measured the cosmic-ray antiproton spectrum and antiproton-to-proton ratio Aguilar et al. (2016a), and are searching for cosmic-ray antimatter nuclei Collaboration (2018). Measurements such as these provide a powerful way to search for new physics, including the annihilation or decay of dark mater particles in the halo of the Milky Way.

Over the past several years, a number of groups have reported the presence of an excess of $\sim$10-20 GeV antiprotons, consistent with the annihilation products of dark matter Hooper et al. (2015); Cirelli et al. (2014); Bringmann et al. (2014); Cuoco et al. (2017a); Cui et al. (2017). Moreover, an excess of GeV-scale gamma-rays from the Inner Galaxy has now been firmly confirmed Hooper and Goodenough (2011); Hooper and Linden (2011); Abazajian and Kaplinghat (2012); Gordon and Macias (2013); Daylan et al. (2016); Calore et al. (2015a); Ajello et al. (2016). Although these two signals are sensitive to very different systematic uncertainties and backgrounds, it is intriguing that they could both be explained by a $\sim$60 GeV dark matter particle with an annihilation cross section near that predicted for a generic thermal relic Daylan et al. (2016); Calore et al. (2015b); Agrawal et al. (2015); Achterberg et al. (2015); Cuoco et al. (2017b). The concordance between these two signals is suggestive and provides considerable motivation for additional indirect searches for dark matter (see e.g. Refs. Geringer-Sameth et al. (2015); Albert et al. (2017)).

Taken at face value, the statistical significance of the AMS-02 antiproton excess is quite high. The authors of Ref. Cuoco et al. (2017a), for example, assess the significance of this signal to be $\sim$ 4.5$\sigma$. It is somewhat challenging, however, to quantify the systematic uncertainties associated with this excess, including, i) uncertainties associated with the injected spectra of cosmic-ray protons, helium, and heavier nuclei, ii) uncertainties associated with the propagation of cosmic rays through the interstellar medium (ISM), iii) uncertainties associated with the antiproton production cross sections for proton-proton, proton-nucleon and nucleon-nucleon inelastic collisions, and iv) uncertainties associated with the impact of the solar wind on the cosmic-ray spectra observed at Earth Cholis et al. (2017).

The remainder of this article is structured as follows. In Sec. II we describe our treatment of the systematic uncertainties listed in the previous paragraph. We then present our main results in Sec. III, finding that the antiproton-to-proton ratio measured by AMS-02 favors the presence of a contribution from annihilating dark matter at the level of 4.7$\sigma$. In Sec. IV, we discuss how the stochastic acceleration of antiprotons in supernova remnants can impact our results, favoring a similar range of dark matter models but with a somewhat lower statistical significance of 3.3$\sigma$. The inclusion of this contribution also leads to a better fit to the antiproton spectrum at energies above $\sim$100 GeV. In Sec. V we summarize our results and discuss them within the larger context of indirect searches for annihilating dark matter.

II Methodology and Assumptions

In this section, we describe our efforts to quantify the systematic uncertainties associated with the cosmic-ray antiproton spectrum. Readers interested only in the results of our analysis may feel free to skip to Sec. III.

II.1 Cosmic-Ray Injection and Propagation in

the Interstellar Medium

Antiprotons can be produced as secondary cosmic rays when energetic cosmic-ray primaries (*i.e. *cosmic rays accelerated by supernova remnants) collide with interstellar gas. Cosmic rays that acquire their energy through first order Fermi acceleration are generally expected to exhibit a power-law spectrum, $dN/dE\propto E^{-\alpha}$, with a typical spectral index of $\alpha\sim 2.0$ (see *e.g. *Ref. Mertsch (2018)). A large number of supernova remnants contribute to the cosmic-ray spectrum, with a variety of ages and at a range of distances from the Solar System. In our calculations, we adopt the following parameterization for the average injected spectra of cosmic-ray protons and nuclei:

[TABLE]

where $R$ is the cosmic-ray rigidity (the momentum-to-charge ratio). For simplicity, we adopt the same values of $\alpha_{1}$, $\alpha_{2}$ and $E_{\rm br}$ for protons, helium and other nuclear species, and do not account for the spectral hardening that has been observed at rigidities above $\simeq$ 200 GV in proton, He, Li, Be, B, C, and O cosmic rays Aguilar et al. (2015a, b, 2017, 2018). While ignoring the possibility of a spectral break above 200 GV could impact the antiproton flux by $\simeq 1\%$ at 1 GV (2%, 20% at 10 and 100 GV, respectively), such a feature cannot produce any spectral features in the $\sim$10-20 GeV range and has a negligible impact on our results. Furthermore, given that the B/C ratio shows no evidence of such hardening Aguilar et al. (2016b), we conclude that this high-energy spectral feature is likely to be the result of variations in the local source distribution, not unlike that observed in the spectrum of cosmic-ray positrons produced by pulsars (see e.g. Refs. Cholis et al. (2018a, b)).

Once injected into the ISM, cosmic rays undergo diffusion, convection, and diffusive reacceleration (more massive nuclei may also experience fragmentation). Energy losses for cosmic-ray nuclei by ionization and Coulomb collisions are also included but have a subdominant impact on our results. We take each of these processes into account by solving the cosmic-ray transport equation numerically, using the publicly available code Galprop v54 1.984 v54 availabe at: https://gitlab.mpcdf.mpg.de/aws/galprop ; http://galprop.stanford.edu ; Strong (2015); Strong and Moskalenko (1998).

To model the effects of isotropic and homogenous transport throughout a zone extending up to a half-height of $z_{L}$ from the Galactic Disk, we adopt the following diffusion coefficient:

[TABLE]

where $\beta\equiv v/c$ and $\delta$ is the diffusion index associated to the spectral index of magnetohydrodynamic turbulence in the ISM. Typical values for $\delta$ are 0.33 for Kolmogorov turbulence Kolmogorov (1941) and 0.5 for the Kraichnan case Kraichnan and Montgomery (1980). Values within this range are generally consistent with the existing body of cosmic-ray data (see e.g. Ref. Trotta et al. (2011)).

Diffusive reacceleration is described by a diffusion coefficient in momentum space Seo and Ptuskin (1994),

[TABLE]

where the Alfv$\acute{\textrm{e}}$n speed, $v_{A}$, is the speed at which hydromagnetic waves propagate in the ISM plasma.

The convective wind speed, $v_{c}$, has a gradient perpendicular to the Galactic Plane,

[TABLE]

To constrain the above parameters, we follow the procedure described in Refs. Cholis et al. (2016, 2017), using data from AMS-02, PAMELA and Voyager 1. Previous studies have found the $\sim$10-20 GeV antiproton excess to be robust to variations in the values of these parameters Cuoco et al. (2017a); Cui et al. (2017); Hooper et al. (2015). Instead of repeating the same procedure here, we chose to adopt three representative models for cosmic-ray injection and transport (see Table 1), each of which provides a good overall fit to the observed cosmic-ray proton, helium, carbon, and boron-to-carbon ratio spectra up to 200 GV (see also Table I of Ref. Cholis et al. (2017)).111Our ISM models ”I”, ”II” and ”III” are the same as ISM models ”C”, ”E and ”F” from our previous work in Ref. Cholis et al. (2016, 2017). We find that if we add the ISM models ”A”, ”B” and ”D” from Ref. Cholis et al. (2016), our results regarding the GeV $\bar{p}/p$ excess at the 5-20 GeV in $E_{\textrm{kin}}$ fall within the same range of significance.

II.2 The Antiproton Production Cross Section

The production of antiprotons in the inelastic collisions of high-energy nuclei has been studied at a number of collider experiments Dekkers et al. (1965); Capiluppi et al. (1974); Allaby et al. (1970); Guettler et al. (1976); Johnson et al. (1977); Antreasyan et al. (1979); Arsene et al. (2007); Anticic et al. (2010); Acharya et al. (2018), and this information has been parameterized in several different forms Tan and Ng (1982, 1983); Duperray et al. (2003); di Mauro et al. (2014); Kappl and Winkler (2014); Kachelriess et al. (2015); Bernstein (2018). Most of the existing work on this subject has focused on the direct productions of antiprotons in proton-proton collisions. Also relevant to the problem at hand is the production of antineutrons, which lead to the production of antiprotons through their decay in the ISM. Furthermore, significant uncertainties apply to the rate of antiproton production from helium and other nuclei, which are collectively responsible for approximately 40% of the overall flux.

Using Galprop, we account for the production of antiprotons from all cosmic-ray species Moskalenko et al. (2002). Based on Ref. di Mauro et al. (2014), the 3$\sigma$ uncertainty on $\sigma_{pp\rightarrow X+\bar{p}}$ is about $\simeq 40\%$ for antiprotons with a kinetic energy of 1 GeV. This uncertainty is energy dependent, as shown in Fig. 8 of Ref. di Mauro et al. (2014). Following up on our previous work Cholis et al. (2017), we first evaluate the antiproton flux for a given cosmic-ray transport model and cross section, and then marginalize over a flat prior within the energy-dependent $3\sigma$ uncertainties on the antiproton production cross section, as quoted in Ref. di Mauro et al. (2014). We do this through the following energy-dependent scaling factor:

[TABLE]

The parameters $a$, $b$ and $c$ and are allowed to vary over a large range of values in order to cover the range quoted in Ref. di Mauro et al. (2014). By allowing this energy dependent cross section to vary without penalty within the quoted 3$\sigma$ uncertainties, we are conservatively allowing for a generous range of behavior in our analysis. By including a greater degree of flexibility in this parameterization (*e.g. *a term proportional to $[\ln(E_{\rm kin}/{\rm GeV})]^{3}$), one could absorb much of the $\sim$10-20 GeV excess observed in the cosmic-ray antiproton spectrum. We emphasize, however, that the terms contained in Eq. 5 more than adequately encompass the physically plausible range of uncertainties associated with this quantity. Moving forward, high precision laboratory measurements of the antiproton production cross section could reduce these uncertainties and substantially increase our ability to search for dark matter annihilation products in the cosmic-ray spectrum.

II.3 Solar Modulation

As cosmic rays enter the Solar System, their spectra are modulated by the solar wind and its embedded magnetic field. We use the standard formula to model the impact of the modulation potential Gleeson and Axford (1968):

[TABLE]

where $E_{\textrm{kin}}$ is the kinetic energy of the cosmic ray measured at Earth, Z$e$ and $m$ are the charge and mass of the cosmic ray, $dN^{\oplus}/dE_{\textrm{kin}}$ is the spectrum measured at Earth and $dN^{\rm ISM}/dE_{\rm kin}^{\rm ISM}$ is the spectrum present in the ISM, prior to the effects of solar modulation. Based on Ref. Cholis et al. (2016), we adopt the following rigidity, time and charge-dependent modulation potential:

[TABLE]

where $B_{\rm tot}(t)$ is the strength of the heliospheric magnetic field at Earth (as measured by ACE http://www.srl.caltech.edu/ACE/ASC/ ), $A(t)$ is its polarity, and $\alpha(t)$ is the tilt angle of the heliospheric current sheet (based on models provided by the Wilcox Solar Observatory http://wso.stanford.edu/Tilts.html ). $R$ is the rigidity of the cosmic ray prior to entering the Solar System, and $R_{0}\equiv 0.5$ GV.

To model the uncertainties associated with solar modulation, we allow for $0.32\leq\phi_{0}\leq 0.38$ GV and $0\leq\phi_{1}\leq 16$ GV, each of which represent the 2$\sigma$ range presented in Ref. Cholis et al. (2017). The quantity $N^{\prime}(q)$, along with averaged values of $B_{\rm tot}(t)$ and $\alpha(t)$, are given in Table II of Ref. Cholis et al. (2017) for each six-month interval. We perform the Solar modulation correction over each of the these six-month intervals and calculate the averaged spectra before comparing to the data (see Refs. Cholis et al. (2016, 2017) for further details).

III Results

III.1 Fitting the Antiproton-to-Proton Ratio Without Dark Matter

We begin by considering the antiproton-to-proton ratio presented by the AMS-02 Collaboration Aguilar et al. (2016a), performing the fit without any contribution from dark matter or other exotic physics. We treat the astrophysical and particle physics uncertainties as described in the previous section, and consider each model for the injection and propagation in the ISM independently (see Table 1). For each ISM model, we scan across a six-dimensional grid ($\phi_{0}$, $\phi_{1}$, $a$, $b$, $c$, and the normalization of the ISM gas density), calculating the log-likelihood for each point in this parameter space and then selecting the combination of values that provides the best fit to the data. The results of our fit are shown in Fig. 1. The grey band shown in each frame represents the combined uncertainties associated with solar modulation and the antiproton production cross section. Although these bands are quite wide, we emphasize that these uncertainties are highly correlated, and are not generally capable of producing (or absorbing) narrow spectral features, such as those which might arise from annihilating dark matter. In each case, our model provides a reasonably good description of the data, yielding a $\chi^{2}$ per degree-of-freedom of 0.79, 1.28 and 1.03 for ISM Models I, II and III, respectively (as these error bars include both statistical and systematic uncertainties, we caution that one cannot use these values to formally address the quality of the fit). In the lower frames of this figure, however, one can easily identify a positive residual which appears at $\sim$10-20 GeV (and a deficit at $\sim$5-10 GeV), as well as an excess at energies above $\sim$100 GeV.

III.2 Including Annihilating Dark Matter

The spectrum of antiprotons produced in dark matter annihilation processes can be calculated using Monte Carlo event generators such as PYTHIA Sjostrand et al. (2008) and HERWIG Corcella et al. (2001). In this study, we use the publicly available PPPC4DMID code Cirelli et al. (2011) which provides the differential spectra of antiprotons from DM annihilations, $dN_{\bar{p}}/dE_{\bar{p}}$. Although throughout most of this study we focus on the representative case of annihilations to $b\bar{b}$, we consider in the Appendix models in which the dark matter annihilates to light quarks or to $W^{+}W^{-}$. The PPPC4DMID code includes electroweak corrections which are important in the case of heavy dark matter particles, when the annihilation products can be highly boosted and emit a $W$ or $Z$ before decaying or hadronizing Ciafaloni et al. (2011).

For the distribution of dark matter in the Milky Way we adopt an Navarro-Frenk-White (NFW) profile Navarro et al. (1996):

[TABLE]

We set the normalization parameter, $\rho_{0}$, such that the local density (at $r=8.5$ kpc) is 0.4 GeV/cm3 Catena and Ullio (2010); Salucci et al. (2010) and adopt a scale radius of $r_{s}=20$ kpc. We note that the results presented here are not highly sensitive to the choice of the halo profile. If we had instead adopted an Einasto profile Einasto (1965) or a profile with a slightly steeper inner slope (as motivated by the observed profile of the Galactic Center gamma-ray excess Daylan et al. (2016); Calore et al. (2015a)), the local antiproton spectrum would be largely unaffected. The reason for this is that most of the cosmic rays in the energy range of interest originate from the surrounding few kpc, and thus the dependence on the dark matter halo profile is largely limited to the overall normalization (*i.e. *the local density).

In Fig. 2, we show the impact of annihilating dark matter on the fit to the antiproton-to-proton ratio for the case of dark matter annihilating to $b\bar{b}$. For the case of the thermal relic benchmark cross section (shown as a white dashed line Steigman et al. (2012)), this data excludes (at the $2\sigma$ level) dark matter masses up to 47 GeV and between 136-286 GeV, representing one of the strongest constraints on annihilating dark matter. There are two regions of parameter space, however, in which a dark matter annihilation signal improves the quality of the fit. The best overall fit is found for the case in which a $m_{\chi}=64-88$ GeV dark matter candidate annihilates with a cross section of $\sigma v=(0.8-5.2)\times 10^{-26}$ cm3/s. Such a contribution improves the fit by $2\Delta\ln\mathcal{L}=22.0$, 59.8 and 54.2 for ISM Models I, II, and III, respectively, corresponding to a statistical preference between 4.7 and 7.7$\sigma$. It is noteworthy how similar these parameters are to those that are required to generate the observed characteristics of the Galactic Center gamma-ray excess Daylan et al. (2016); Calore et al. (2015a). At higher masses (\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}1 TeV), annihilating dark matter particles can also improve the fit to this dataset, although to a lesser extent. We remind the reader that at each point in the fit we have marginalized over the parameters associated with the antiproton production cross section and solar modulation as described in Sec. II, and therefore our results indicate that the presence of this excess is statistically significant, even in light of these systematic uncertainties.

In Fig. 3, we show the spectrum of the antiproton-to-proton ratio, including the best-fit contribution from annihilating dark matter. The residual plots (lower frames) clearly illustrate the preference for a contribution from annihilating dark matter peaking in at energies near $\sim$10-20 GeV. In the top three rows of Table 2 we summarize our results, listing the values of the dark matter mass and annihilation cross section that are favored by this fit, for each of the three cosmic-ray injection and transport models considered in this study. In each case, we find a statistically significant preference for a contribution from annihilating dark matter.

We note that our analysis arrives at qualitatively different conclusions than those presented in Ref. Reinert and Winkler (2018), which finds that the statistical significance of the antiproton excess can be reduced to approximately 2.2$\sigma$ after systematic uncertainties are taken into account. We note that there are many significant differences between the cosmic-ray propagation models employed between these papers. Most notably, the authors of Ref. Reinert and Winkler (2018) utilize an analytic two-zone cosmic-ray propagation model, with parameters that are tuned to the antiproton data, as well as to the cosmic-ray positron flux. We utilize numerical cosmic-ray propagation models based on the Galprop code, and choose not to normalize our astrophysical background models to cosmic-ray antiprotons (to avoid biasing the results) or to the spectrum of cosmic-ray leptons (which have vastly different cooling times). In addition, there are significant differences in our modeling of the antiproton production cross section and in our treatment of solar modulation.

Up to this point, we have not considered the possibility that secondary antiprotons could be accelerated in the environments surrounding supernova remnants Cholis et al. (2017); Blasi (2009); Mertsch and Sarkar (2009); Cholis and Hooper (2014); Kohri et al. (2016); Fujita et al. (2009). In the following section we will consider how such a contribution could impact our results.

IV Stochastic Acceleration of Secondary Cosmic Rays in Supernova Remnants

In the standard picture, cosmic rays are produced when a supernova shockfront expands and sweeps through the ISM, trapping particles within its turbulent magnetic field structure long enough for them to be accelerated. These particle species can also interact with the dense gas on either side of the shockfront. Cosmic rays undergo inelastic scattering and decay at the following rate:

[TABLE]

where $\sigma^{\rm inelastic}_{i}$ and $\tau^{\rm dec}_{i}$ are the inelastic scattering cross section and lifetime of cosmic ray species, $i$, and $n_{\rm gas}$ is the number density of gas. If the timescale for acceleration is much shorter than that of inelastic scattering or decay, the secondaries will be efficiently accelerated. Following Refs. Blasi (2009); Mertsch and Sarkar (2009); Cholis and Hooper (2014); Ahlers et al. (2009), we assume Bohm diffusion for the cosmic rays near the shockfront:

[TABLE]

where $r_{L}$ the Larmor radius of the cosmic rays within the magnetic fields and $K_{B}\simeq(B/\delta B)^{2}$ Blasi (2009) quantifies the turbulent nature of the magnetic fields around the shockfront.

The contribution to the cosmic-ray antiproton spectrum from secondary acceleration depends on the value of $K_{B}$ as well as the density of gas in the scattering region, $n_{\rm gas}$. Both of these parameters have a similar impact on the resulting antiproton spectrum, with larger values leading to a higher antiproton-to-proton ratio at high energies. Increasing $K_{B}$ or $n_{\rm gas}$ will also increase the boron-to-carbon ratio at high energies, and this information can be used to independently constrain the values of these parameters Cholis and Hooper (2014); Cholis et al. (2017); Tomassetti and Oliva (2017).

We begin by adopting values for these parameters that provide a good fit to the observed boron-to-carbon ratio: $n_{\rm gas}=2.0$ cm*-3* and $K_{B}=3.05$, 5.2 and 3.7 for ISM models I, II and III, respectively Cholis and Hooper (2014); Tomassetti and Oliva (2017). The results for these cases are shown in Figs. 4 and 5. The presence of the contribution from accelerated secondaries almost entirely removes the excess at energies above $\sim$100 GeV, erasing the preference for $\sim$1-3 TeV dark matter seen in Fig. 2 Cholis et al. (2017). The evidence for a lighter dark matter particle persists in the presence of accelerated secondaries, however, favoring $m_{\chi}=46-94$ GeV and $\sigma v=(0.7-3.3)\times 10^{-26}$ cm3/s with a statistical significance of 3.3$\sigma$ (see Table 2).

By fixing the values of $K_{B}$ and $n_{\rm gas}$ in our calculations to those which reproduce the observed boron-to-carbon ratio, we are implicitly making the assumption that carbon and protons are accelerated in the same supernova remnants. It is conceivable that carbon nuclei and protons are preferentially accelerated in different subsets of the supernova remnant population, with different average values of $K_{B}$ and $n_{\rm gas}$. With this possibility in mind, we repeat the above fit, allowing the impact of secondary acceleration to vary. In particular, we set $K_{B}=6.1$, 10.4 and 7.4 for ISM models I, II and III, respectively, and allow the value of $n_{\rm gas}$ to float freely in the fit (see Table IV of Ref. Cholis et al. (2017)). We show the results of this fit in Figs. 6 and 7. Although the value of $n_{\rm gas}$ takes on different values throughout the mass-cross section plane, we note that our best-fit points correspond to $n_{\rm gas}=0.39$ cm*-3*, 1.28 cm*-3* and 1.57 cm*-3* for ISM models I, II and III, respectively. The fact that these values are similar to those favored to explain the measured boron-to-carbon ratio suggests that carbon nuclei and protons are likely accelerated in the same class of astrophysical sources.

Somewhat surprisingly, the favored parameter space in the right frame of Fig. 6 is entirely surrounded by an excluded region. This stems from the competing demands of the fitting algorithm to match both the (1) high-energy (\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}100 GeV) antiproton excess, which requires large values of $n_{\rm gas}$, and the (2) low-energy ($\sim$10-20 GeV) antiproton excess, which is best fit by dark matter in a scenario with a smaller value of $n_{\rm gas}$. For dark matter masses near our best fit value, moderate values of the annihilation cross section are disfavored, because they force the fit to overproduce the low-energy antiproton excess for the values of $n_{\rm gas}$ that provide the best fit to the high-energy data.

V Discussion and Summary

In this article, we have studied the cosmic ray antiproton-to-proton ratio measured by AMS-02 Aguilar et al. (2016a), and considered the implications of this measurement for dark matter annihilating in the halo of the Milky Way. Our main results are summarized in Table 2 and in Fig. 8. In each case considered, we find a significant excess of $\sim$10-20 GeV antiprotons, even after marginalizing over a generous range of parameters associated with the effects of solar modulation and the antiproton production cross section. This excess is well fit by annihilating dark matter particles, with a mass and cross section in the range of $m_{\chi}\approx 46-94$ GeV and $\sigma v\approx(0.7-5.2)\times 10^{-26}$ cm${}^{3}/$s, respectively (for the representative case of annihilations to $b\bar{b}$). Other annihilation channels can also provide a good fit, although for slightly different parameter ranges (see Appendix).

Although this result is interesting in its own right, it is particularly intriguing that the range of dark matter models that can accommodate the antiproton excess is very similar to those which could generate the excess of GeV-scale gamma rays observed from the Galactic Center Hooper and Goodenough (2011); Hooper and Linden (2011); Abazajian and Kaplinghat (2012); Gordon and Macias (2013); Daylan et al. (2016); Calore et al. (2015a); Ajello et al. (2016). In the left frame of Fig. 8 we compare the regions of dark matter parameter space that are able to account for the gamma-ray excess Calore et al. (2015a) to those favored by the analysis of the antiproton spectrum presented in this study. These two regions overlap, and collectively favor dark matter particles with $m_{\chi}=48-67$ GeV and $\sigma v=(1.4-2.4)\times 10^{-26}$ cm${}^{3}/$s.

Putting the antiproton excess aside for a moment, our analysis also yields stringent constraints on the dark matter annihilation cross section, in many cases competitive with, or more stringent than, other bounds. In the right frame of Fig. 8, we show our overall constraint on the dark matter annihilation cross section, which we take to be the weakest of the constraints shown in Figs. 2, 4 and 6, evaluated at each value of $m_{\chi}$. Compared to the constraints derived from gamma-ray observations of dwarf spheroidal galaxies Albert et al. (2017), we find that the limit presented in this study is stronger for dark matter particles with a mass below 40 GeV or between 130 and 540 GeV (for annihilations to $b\bar{b}$).

As this article was being finalized, Ref. Cuoco et al. (2019) appeared on the arXiv which addresses many of the same questions discussed here. The authors of Ref. Cuoco et al. (2019) reach conclusions that are very similar to our own.

Acknowledgements.

We would like to thank Simeon Bird and Marc Kamionkowski for valuable discussions. DH is supported by the US Department of Energy under contract DE-FG02-13ER41958. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE- AC02-07CH11359 with the US Department of Energy.

Appendix A Results For Other Dark Matter Annihilation Channels

Throughout this study, we have focused on the representative case of dark matter particles that annihilate to $b\bar{b}$. Clearly this is not the only possibility, and dark matter that annihilates to other final states could also be responsible for the antiproton excess observed by AMS-02. In Fig. 9 and Table 3 we show our results for dark matter candidates that annihilate to light quarks or to $W^{+}W^{-}$. We also note that hidden sector dark matter candidates could produce a similar spectrum of antiprotons, in particular within the context of Higgs portal models (see Fig. 12 of Ref. Escudero et al. (2017)). In such a model, the dark matter would annihilate to other hidden sector states, which then decay through mixing with the Standard Model Higgs boson. We note that we utilized PYTHIA Sjostrand et al. (2008) for dark matter particles lighter than 86 GeV in the $W^{+}W^{-}$ case, as PPPC Cirelli et al. (2011) does not generate reliable results in the case of $m_{\chi}\approx m_{W}$.