Investigating the dark matter signal in the cosmic ray antiproton flux with the machine learning method

TL;DR

This study uses machine learning to analyze cosmic ray antiproton data, revealing that the detection of dark matter signals strongly depends on the interaction and propagation models, with some models indicating potential signals around 100 GeV to 1 TeV.

Contribution

It introduces a machine learning approach to efficiently explore complex parameter spaces in dark matter signal analysis from cosmic ray data.

Findings

Most models show no significant dark matter signal above 2σ.

The EPOS-LHC model suggests a >3σ signal at ~1 TeV.

A significant signal at ~100 GeV is found in the diffusive reacceleration model, but diminishes with charge-dependent solar modulation.

Abstract

We investigate the implications on the dark matter (DM) signal from the AMS-02 cosmic antiproton flux. Global fits to the data are performed under different propagation and hadronic interaction models. The uncertainties from the injection spectrum, propagation effects and solar modulation of the cosmic rays are taken into account comprehensively. Since we need to investigate extended parameter regions with multiple free parameters in the fit, the machine learning method is adopted to maintain a realistic time cost. We find all the effects considered in the fitting process interplay with each other, among which the hadronic interaction model is the most important factor affecting the result. In most hadronic interaction and CR propagation models no DM signal is found with significance larger than $2\sigma$ except that the EPOS-LHC interaction model requires a more than $3\sigma$ DM signal with DM mass around $1\,\mathrm{TeV}$. For the diffusive reacceleration propagation model there is a highly significant DM signal with mass around $100\,\mathrm{GeV}$. However, the signal becomes less than $1\sigma$ if we take a charge dependent solar modulation potential in the analysis.

Tables

Table 1. Table 1: The best fits for the propagation parameters and the positive modulation potential ϕ + subscript italic-ϕ \phi_{+} in the DR and DR-2 models, obtained from the Ref.  Yuan  et al. ( 2017 ) .

	$D_0$	$\delta$	$z_h$	$v_A$	$\eta$		${\varphi }_{+}$
	${10}^{28}{\mathrm{cm}}^2{s}^{-1}$		$\mathrm{kpc}$	$\mathrm{km}{s}^{-1}$			$\mathrm{GV}$
DR	7.24	0.38	5.93	38.5	1		0.527
DR-2	4.16	0.5	5.02	18.4	-1.28		0.636

Equations

$$\frac{\partial \psi }{\partial t}$$

$$f(r,z)={\left(\frac{r}{r_{\odot }}\right)}^a\exp \left[-\frac{b(r-r_{\odot })}{r_{\odot }}\right]\exp \left(-\frac{\left|z\right|}{z_s}\right),$$

$$D_{xx}=D_0{\beta }^{\eta }{\left(R/R_0\right)}^{\delta }.$$

$$\begin{gathered} D_{pp}{D}_{xx}=\frac{4p^2{v}_A^2}{3\delta (4-{\delta }^2)(4-\delta )\omega } \\ , \end{gathered}$$

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

$$\inf \left\{\frac{1}{n}\sum _{i=1}^n\max (0,1-y_i(\mathbf{\omega }\cdot {\mathbf{x}}_i-b)+\lambda |\mathbf{\omega }{|}^2|\mathbf{\omega }\in {\mathcal{R}}^l,b\in \mathcal{R}\right\},$$

$$\inf \left\{U\sum _{i=1}^n\max \left(\epsilon ,|\mathbf{\omega }\cdot {\mathbf{x}}_i+b-y_i|\right)+\frac{1}{2}|\mathbf{\omega }{|}^2|\mathbf{\omega }\in {\mathcal{R}}^l,b\in \mathcal{R}\right\},$$

$$\hat{y}(\phi )=\mathbf{W}\cdot \mathbf{\phi }+b(\mathbf{\phi }\in {\mathcal{R}}^k,k>l).$$

$$k({\mathbf{x}}_i,{\mathbf{x}}_j)\equiv \hat{\mathbf{\phi }}({\mathbf{x}}_i)\cdot \hat{\mathbf{\phi }}({\mathbf{x}}_j).$$

$$k({\mathbf{x}}_i,{\mathbf{x}}_j)=\exp \left(-\gamma |{\mathbf{x}}_i-{\mathbf{x}}_j{|}^2\right),$$

$$\inf \left\{U\sum _{i=1}^n\max \left(\epsilon ,|k(\mathbf{\omega },{\mathbf{x}}_i)+b-y_i|\right)+\frac{1}{2}k(\mathbf{\omega },\mathbf{\omega })|\mathbf{\omega }\in {\mathcal{R}}^l,b\in \mathcal{R}\right\}.$$

$$\begin{array}{rlr}\text{minimize}& & U\sum _{i=1}^N{\xi }_i^{\ast }+U\sum _{i=1}^N{\xi }_i+\frac{1}{2}k(\mathbf{\omega },\mathbf{\omega })\\ \text{subject to}& & \{\begin{array}{rl}{y}_i-k(\mathbf{\omega },{\mathbf{x}}_i)-b\le & \epsilon +{\xi }_i^{\ast }\\ -y_i+k(\mathbf{\omega },{\mathbf{x}}_i)+b\le & \epsilon +{\xi }_i\\ {\xi }_i^{\ast }\ge 0{\xi }_i\ge 0& \end{array}\end{array}.$$

$$\mathcal{L}=$$

$$-\sum _{i=1}^N{\alpha }_i(y_i-\mathbf{W}\cdot \hat{\mathbf{\phi }}({\mathbf{x}}_i)-b+\epsilon +{\xi }_i)-\sum _{i=1}^N{\gamma }_i({\xi }_i^{\ast }+{\xi }_i),$$

$$\begin{array}{rl}\text{minimize }& \epsilon \sum _{i=1}^n({\alpha }_i^{\ast }+{\alpha }_i)+\sum _{i=1}^n{y}_i({\alpha }_i-{\alpha }_i^{\ast })-\frac{1}{2}\sum _{i,j=1}^N({\alpha }_i-{\alpha }_i^{\ast })k({\mathbf{x}}_i,{\mathbf{x}}_j)({\alpha }_j-{\alpha }_j^{\ast })\\ \text{subject to }& \{\begin{array}{rl}& {\alpha }_i={\alpha }_i^{\ast }\\ & 0\le {\alpha }_i,{\alpha }_i^{\ast }\le U\end{array}\end{array}.$$

Full text

I Introduction

Cosmic ray (CR) antiproton is a very promising probe in the dark matter (DM) indirect detection. In recent years, a great progress of the measurement has been made by the AMS-02 experiment, which precisely measured the CR $\bar{p}/p$ ratio and $\bar{p}$ flux up to $\sim 450$ GeV. Based on these results, lots of theoretical studies have been performed to investigate the possible contribution from DM in the antiproton flux Giesen et al. (2015); Jin et al. (2015); Chen et al. (2015); Ibe et al. (2015); Hamaguchi et al. (2015); Kohri et al. (2015); Kappl et al. (2015); Lu and Zong (2015); Lin et al. (2015); Cui et al. (2016); Cuoco et al. (2017); Feng et al. (2016); Huang et al. (2017); Lin et al. (2016); Cui et al. (2018a, b); Cuoco et al. (2019); Cholis et al. (2019). As the deviation between the expected secondary CR antiproton flux produced by the astrophysical processes and the data is not significant, any slight change of the background expectation could influence the DM implication. Therefore, the theoretical uncertainties of the secondary CR antiproton flux, mainly arising from the CR propagation process and the hadronic interactions between the CR particles and interstellar medium (ISM), are very important in the determination of the DM contribution.

Although the principle of the strong interaction has been well described by the quantum chromodynamics (QCD), only the processes with large momentum transfers can be perturbatively calculated. The production processes of secondary CR antiprotons, which involve the forward scattering processes with multi-particle production, are always calculated under extra simplified assumptions or using some empirical parametrizations. This kind of processes has been a long time difficulty in the cosmic ray study. Many phenomenological models with different simplified assumptions have been constructed in the literature Kalmykov et al. (1997); Ostapchenko (2006a); Engel et al. (1992); Engel (1995); Bopp et al. (2008); Werner et al. (2006); Pierog et al. (2015). In addition, several empirical parametrizations only aiming to reproduce the antiproton production cross section have also been developed Tan and Ng (1983); di Mauro et al. (2014); Winkler (2017). On the other hand, the CR propagation can be determined by the fit to the observed secondary-to-primary nuclei ratios like B/C and (Sc+Ti+V)/Fe, and unstable-to-stable ratios of secondary nuclei like ${}^{10}\mathrm{Be}/^{9}\mathrm{Be}$ and ${}^{26}\mathrm{Al}/^{27}\mathrm{Al}$ Strong and Moskalenko (1998); Gaisser (1990); Di Bernardo et al. (2010); Maurin et al. (2001). With the constraints of current data, there remains a considerable acceptable region for the varying propagation parameters Yuan et al. (2017).

With these uncertainties involved, we have found that several hadronic models can well fit the AMS-02 antiproton data, and derived the upper-limits on the DM annihilation cross section in a previous work Lin et al. (2016). However, in the scenarios with the other hadronic models or astrophysical configurations, the secondary CR antiproton could be insufficient in some energy ranges. In these cases, it is interesting to explore the favored region in the DM parameter space by investigating how an extra DM contribution could improve the fit Cui et al. (2016); Cuoco et al. (2017); Cui et al. (2018a); Cuoco et al. (2019); Cholis et al. (2019).

To set the upper-limits, we could just separately estimate the antiproton fluxes with thousands of benchmark astrophysical parameter sets within the 2$\sigma$ confidence region accounting for the B/C fit, and use the band spanned by the corresponding limits to indicate the astrophysical uncertainties. However, claiming an excess and investigating a favored DM parameter region should be prudent and much more complicated. As the favored region is sensitive to the secondary CR antiproton spectrum, investigating thousands of benchmark astrophysical parameter sets is far from enough. It is necessary to vary the astrophysical parameters when searching for the favored region in order to reach a comprehensive solution Cui et al. (2016); Cuoco et al. (2017); Cui et al. (2018b); Cuoco et al. (2019).

In this work, we try to investigate this favored region by performing a likelihood optimization in the $m_{\chi}-\langle\sigma v\rangle$ plane, in which the profile likelihood for each given $(m_{\chi},\langle\sigma v\rangle)$ is obtained by varying the astrophysical parameters about ten thousand times. Since the CR propagation equation solution for each set of astrophysical parameter always costs about several minutes by using the well-developed numerical packages like GALPROP Strong and Moskalenko (1998); Moskalenko and Strong (1998) or DRAGON Evoli et al. (2017), it is very time-consuming to perform the whole scanning analysis. In order to solve this efficiency problem, we adopt the machine learning method in the analysis, which works through a training and predicting process. We perform the CR propagation calculation using GALPROP for about ten thousand sets of the astrophysical parameters, and take these results to train the machine learning model. Then we utilize the well-trained model to predict the antiproton flux for given astrophysical parameters in the scanning analysis. Since the calculation of the machine learning model is much faster than that of GALPROP, this method is able to significantly improve the efficiency. With this technique introduced, it then is possible to statistically study the DM implication in the AMS-02 antiproton data for different hadronic and CR propagation models.

This paper is organized as follows. In Sec. II, we describe the physical processes related to the CR antiproton flux, including the CR propagation and hadronic models, and DM annihilation. In Sec. III, we introduce the adopted machine learning method in brief, and show the training quality of the machine learning model. In Sec. IV, we show the final scan results for different hadronic and CR propagation models, and discuss the implications of DM contribution in the AMS-02 antiproton flux. Finally, we give the summary in Sec. V.

II Cosmic ray antiproton flux

The CR antiprotons are considered to originate from the hadronic interaction between primary CR particles and ISM, perhaps along with the additional contribution from DM annihilation/decay. These antiprotons would undergo the processes like propagation and solar modulation before arriving at the Earth.

In this section, we would introduce the CR propagation, solar modulation, hadronic interaction, and DM annihilation processes considered in this work in detail. The relevant adopted parameters of the CR propagation and solar modulation processes are based on our previous work Yuan et al. (2017).

II.1 Propagation of Galactic cosmic rays

The Galactic CRs are believed to originate from the supernova remnants (SNRs), and propagate diffusively in the Galactic magnetic field (GMF) that is extended within a cylindrical propagation halo. Such a diffusive propagation together with the accompanying effects can be described with the equation Strong et al. (2007)

[TABLE]

where $\psi$ is the CR density per momentum interval, $Q(\mathbf{x},p)$ is the CR source term which is assumed to follow a three-pieces broken power-law with respect to the rigidity, the $D_{xx}$ term describes the diffusion effect, the $\mathbf{V}_{c}$ terms describe the convection effect, the $D_{pp}$ term describes the reacceleration effect, the term with $\dot{p}\equiv\mathrm{d}p/\mathrm{d}t$ is the momentum loss term, and the time scales $\tau_{f}$ and $\tau_{r}$ characterize fragmentation processes and radioactive decays, respectively.

This equation applies to all kinds of CR particles. As the source term for the secondary particles depends on the fragmentation and decay of the primary particles, the secondary fluxes linearly depend on the primary fluxes. Therefore, their ratios could be used to characterize the propagation process. In practice, the secondary-to-primary ratios like B/C and the unstable-to-stable ratios like ${}^{10}\mathrm{Be}/^{9}\mathrm{Be}$ are usually adopted to constraint the coefficients included in Eq. 1 Putze et al. (2010); Trotta et al. (2011). In the previous work Yuan et al. (2017), we have fitted the B/C ratio observed by AMS-02 Aguilar et al. (2016a) with the different propagation models, and found that the diffusion plus reacceleration (DR) model works better than the diffusion plus convection (DC) model. Therefore, we only consider the reacceleration effect in this work.

The spatial distribution of the primary CR sources $Q(\mathbf{x},p)$ is supposed to follow the SNR distribution

[TABLE]

where $r_{\odot}=8.5\,\mathrm{kpc}$ and $z_{s}\simeq 0.2\,\mathrm{kpc}$ are the distance from the Sun to the Galactic Center and the characteristic height of the Galactic disk, respectively. In this work, we set $a=1.25$ and $b=3.56$ based on the Fermi-LAT gamma-ray data following the Ref. Trotta et al. (2011).

The diffusion coefficient $D_{xx}$ of the CR particles is empirically parameterized as a power-law function of rigidity $R$ along with a low energy modification Maurin et al. (2010)

[TABLE]

The power-index $\delta$ is expected to be $1/3$ for a Kolmogorov spectrum of interstellar turbulence, or $1/2$ for a Kraichnan cascade. However, we just treat it as a free parameter and determine it by the observed B/C ratio. $\beta$ is the CR particle velocity in units of the light speed that would tends to 1 in the relativistic limit. The factor $\beta^{\eta}$ is introduced to describe the effect of turbulence dissipation, where $\eta$ is assumed to be 1 by default. But it is found that changing the value of $\eta$ could help to improve the fitting of secondary CR particles Di Bernardo et al. (2011). Therefore, in this work we also discuss the propagation model with a free $\eta$, which is referred to as the DR-2 model following the Ref. Yuan et al. (2017).

The CR reacceleration effect is a second-order Fermi acceleration effect caused by the collisions between CR particles and the interstellar random weak hydrodynamic waves. It is described by a diffusion in the momentum space whose diffusion coefficient $D_{pp}$ is related to $D_{xx}$ by Seo and Ptuskin (1994)

[TABLE]

where $v_{A}$ is the Alfvén velocity, and $\omega$ is the energy density ratio between the magneto-hydrodynamic wave and the magnetic field. In practical calculation, $\omega$ could be absorbed into the $v_{A}$.

In all, the relevant propagation parameters include the normalizing factor $D_{0}$, power-index $\delta$, reference rigidity $R_{0}$ for the diffusion coefficient, the Alfvén velocity $v_{A}$, the half height of propagation halo $z_{h}$, and perhaps the parameter $\eta$. In this work, we adopt the DR and DR-2 models from Ref. Yuan et al. (2017), where the reference rigidity $R_{0}$ is taken to be $4\,\mathrm{GV}$. The corresponding best fitting parameters are listed in Table. 1.

II.2 Solar modulation

The CR particles would also be affected by the heliospheric magnetic field (HMF) when they enter the heliosphere. Such an effect that mainly influences the particles with rigidity $R\lesssim 20\,\mathrm{GV}$ is called solar modulation, and can be described by the Parker’s transport equation Parker (1965). This equation was firstly solved by Gleeson et. al. Gleeson and Axford (1968) using the force field approximation (FFA), in which the solar modulation potential $\phi$ is the sole parameter.

However, this solution is derived under the spherical symmetric assumption, and has neglected the drift effect caused by the configuration of HMF. A recent study, which solved the Parker’s transport equation with the realistic simulation Potgieter (2014), has found that this drift effect would indeed lead to a charge-sign dependent behaviour in the CR spectra. Therefore, adopting the FFA with only one modulation potential $\phi$ is insufficient to describe all the CR particles. In this work, we adopt the FFA to deal with the solar modulation, but assume two modulation potentials $\phi_{+}$ and $\phi_{-}$ for the positive and negative charged particles respectively.

The strength of solar modulation would vary with the solar activity, thus the value of modulation potential depends on the observation period of the experiment. In Ref. Yuan et al. (2017), we assumed that the modulation potential $\phi$ is linearly correlated with the sunspots’ number $N$, such that $\phi=\phi_{0}+\phi_{1}\cdot N/N_{\mathrm{max}}$, and fitted the data from different periods for the coefficients $\phi_{0}$, $\phi_{1}$. For the period of AMS-02 $\bar{p}/p$ observation, the corresponding potentials are $\sim 0.527\,\mathrm{GV}$ and $0.636\,\mathrm{GV}$ for the DR and DR-2 models, respectively, as shown in Table. 1. As these values are derived by fitting the proton and B/C data, they are in fact the modulation potential for positive particles $\phi_{+}$. The negative one $\phi_{-}$ is set to be a free parameter vary from $0.5\phi_{+}$ to $1.5\phi_{+}$ in our analysis.

II.3 Hadronic models

The hadronic interaction processes between the primary CR particles and the environment gas could not be directly calculated in the first principle with QCD. Basing on collider data, many phenomenological and empirical models for these processes have been developed to deal with this problem.

The phenomenological models describe the fragmentation process of nucleon based on quantum field theory along with some effective assumption. These models are derived for different purposes. HERWIGCorcella et al. (2001), PYTHIASjostrand et al. (2006) and SHERPAGleisberg et al. (2009) focus on the hard-scattering processes for high energy physics studies. QGSJET01Kalmykov et al. (1997), QGSJET IIOstapchenko (2006a, b) and SIBYLLEngel et al. (1992); Fletcher et al. (1994); Ahn et al. (2009) focus on the bulk production of soft particles for simulating extensive air showers caused by high energy CR particles. PHOJETEngel (1995); Engel and Ranft (1996); Engel et al. (1995), DPMJETBopp et al. (2008) and EPOSWerner et al. (2006); Pierog et al. (2015) fall in between the two categories.

As we aim to calculate the secondary CR particles which are dominantly contributed by the soft production, the models focusing on the hard-scattering processes are not taken into account. The efficiency of the other models in the concerned energy region could be checked by comparing them with the collider data at $\sqrt{s}\sim\mathcal{O}(10\,\mathrm{GeV})$. These data are only available at several benchmark colliding energies in the experiments. The CERN Intersecting Storage Rings (ISR) experiment has measured the $pp\rightarrow\bar{p}X$ cross section at $\sqrt{s}=53\,\mathrm{GeV}$ in mid-1970s Albrow et al. (1973). In recent years, the fixed-target experiment NA49 at CERN Super Proton Synchrotron (SPS) has also measured such a cross section at beam momentum of $158\,\mathrm{GeV}/c$, equaling to $\sqrt{s}=17.27\,\mathrm{GeV}$ in the center-of-mass system Anticic et al. (2010).

In the previous work, we have taken the EPOS, QGSJET and SIBYLL into account and discussed their efficiencies Lin et al. (2016). We approached these models through the common interface provided by the CRMC package Pierog and Ulrich . Both the two versions of EPOS, namely EPOS 1.99 Pierog and Werner (2009) and EPOS LHC Pierog et al. (2015), have been tuned to fit the results from the fixed-target experiments Anticic et al. (2010); Laszlo et al. (2007), and could reasonably reproduce all the measured antiproton cross sections mentioned above. The QGSJETII has only considered the antiproton cross sections in a modified version, namely QGSJETII-04m Kachelriess et al. (2015). This modified model fits the low energy antiproton cross section measurements well, but would always tend to overestimate the low energy secondary CR antiproton flux. The SIBYLL can also fit the antiproton cross sections well, but cannot extend to the extremely low energy region with a beam momentum $\lesssim 60\,\mathrm{GeV}/c$. In all, the only two phenomenological models discussed in this work are chosen to be EPOS 1.99 and EPOS LHC.

Another approach to describe the antiproton production cross sections is the empirical model that fits the measurements with certain analytical parametrization Tan and Ng (1983); Duperray et al. (2003); di Mauro et al. (2014); Winkler (2017). In the following, we select several empirical models as the typical cases, including the default hadronic model embedded in GALPROP by Tan & Ng Tan and Ng (1983), the model fitting to the latest NA49 Anticic et al. (2010) and BRAHMS Arsene et al. (2007) data by di Mauro et al. di Mauro et al. (2014), and the model carefully considering the effects of anti-hyperons and isospin by Winkler Winkler (2017). There are in fact two similar parametrizations in Ref. di Mauro et al. (2014), but we just choose the expression in its Eq. 13 here.

These empirical models only describe the antiproton productions from $pp$ collisions. In order to estimate the secondary CR antiprotons, the contributions of Helium which occupy about 10% of both the CR and the ISM should be considered. Following the procedure of GALPROP GAL (2012), we adopt the code CROSEC by Barashenkov & Polanski to derive the antiproton production cross section of $p$-He and He-He collisions from that of $pp$ collisions.

Although empirical models could always reproduce the accelerator data better compared with phenomenological models at the available collision energies, phenomenological models that constructed with certain physical principles could guarantee a more reasonable estimation when extended to the unexplored collision energy regions. Thus we involve both the two kinds of models in this analysis, including EPOS 1.99, EPOS LHC, and the empirical models by Tan & Ng, di Mauro et al., and Winker.

II.4 Dark matter annihilation

The DM particles could directly annihilate into standard model particles, and produce stable final states like protons/antiprotons, electrons/positrons and photons through the hadronization and decay. Since all the DM annihilation channels to light quarks would lead to similar antiproton spectra, we choose the $b\bar{b}$ channel as a benchmark case in this work. The initial injected antiproton spectra are adopted from PPPC 4 DM ID Cirelli et al. (2011).

The contribution from DM annihilation also depends on the density profile of the Galactic DM halo. Here we adopt the Navarro-Frenk-White (NFW) profile Navarro et al. (1997)

[TABLE]

with $r_{s}=20\,\mathrm{kpc}$ and $\rho_{s}=0.35\,\mathrm{GeV}$. This profile corresponds to a local DM density of $0.4\,\mathrm{GeV}\,\mathrm{cm}^{-3}$ near the Solar system, which is consistent with the recent constraints from the Galactic rotation curve Nesti and Salucci (2013); Benito et al. (2019); Karukes et al. (2019).

III Machine learning

In this work, we attempt to investigate the DM implication through a profile likelihood scan in the $m_{\chi}-\langle\sigma v\rangle$ plane. The profile likelihood for each $(m_{\chi},\langle\sigma v\rangle)$ used in the scan is obtained by marginalizing over the propagation parameters near the best fit shown in Table. 1, which would require us to solve the propagation equation $\sim 10^{4}$ times. Using the numerical package GALPROP Strong and Moskalenko (1998), each propagation solution would cost about one or several minutes. Such a cost is unacceptable for the parameter scan mentioned above, thus a more efficient method is required. We adopt the machine learning method for this purpose.

The supervised machine learning is always used to perform general classifications or regressions. Such regressions could be applied to predict the wanted CR fluxes in order to release us from the time-consuming numerical solutions. We have tested several widely used learning algorithms for this problem, including the support-vector machine (SVM) Cortes (1995), K-nearest neighbors (KNN) Altman (1992), decision tree (DT) Kamiński (2018), random forest Ho (1998), and artificial neural networks (ANN) McCulloch (1943). We find that, focusing on the 2$\sigma$ confidence region favored by the B/C data, the relatively simple SVM algorithm works the best. The reason is that this interested problem is almost a linear regression problem, as the solutions of Eq. 1 could be perturbatively expanded within the narrow region near the best-fit propagation parameters. Therefore, in our scanning analysis, we choose the support-vector regression (SVR) algorithm that based on the SVM. The reader is referred to the Appendix A for the detailed description of this algorithm. We use the package LIBSVM Chang and Lin (2011) to approach training and predicting process of the SVR.

Note that the SVR algorithm can be only applied to the single target problem. However, the propagation solution is described as a multi target problem, whose targets are the CR fluxes at different energies. In order to involve the relation between different targets, we adopt the stacked single-target (SST) algorithm in the calculation Spyromitros-Xioufis (2016). In the SST procedure, we train each target separately, and then combine the predictions of these targets together with the original input to perform a second training. This second training would correct the prediction of each single target with the information of the other targets.

III.1 2$\sigma$ confidence region for the astrophysical parameters

In this work, our scans are constrained within the 2$\sigma$ confidence region derived from the fit to the $\mathrm{B}/\mathrm{C}$, $\mathrm{Be}^{10}/\mathrm{Be}^{9}$, and proton flux observation Yuan et al. (2017). For each propagation model, we adopt a set of samples generated by the Monte Carlo Markov-Chain (MCMC) algorithm from our previous work Ref. Yuan et al. (2017). In each sample, there are a set of propagation, injection, and solar modulation parameters together with the corresponding $\chi^{2}$ for the $\mathrm{B}/\mathrm{C}$ fitting. The best $95\%$ of the samples indicate the 2$\sigma$ confidence region.

In the implementation, such a 2$\sigma$ confidence region is also recognized with machine learning. We perform a training on $\bm{\theta}_{i}\mapsto\chi^{2}_{i}$ with the given samples, and then use the predicting function $\hat{h}_{\chi^{2}}(\bm{\theta})$ to check whether the point $\bm{\theta}$ falls inside the region.

Note that most of the MCMC samples would concentrate around the best fit parameter point and is sparse in the other region. We thus append a set of evenly distributed samples in order to obtain a better prediction in the whole region. In all, there are $20,000\sim 30,000$ MCMC samples and 5000 appended samples for each propagation model. Such combined sample sets would be used as the training sets for all the following analyses.

III.2 Training for the background spectrum

Selecting a sequence of energies which cover the whole energy range of AMS-02, we take the corresponding proton, secondary antiproton and tertiary antiproton fluxes as the training targets. In order to keep the problem simple, we involve no solar modulation effect in the training. As mentioned above, we adopt the SST procedure in which the neighbor fluxes of each point are involved in the training. In Fig. 1, we show a random prediction example for the model DR-2 + EPOS LHC.

It can be seen that the calculation error introduced by this method is obviously smaller than the measurement error, which shows that the introduction of the SVR is sensible in this analysis.

In order to further guarantee the validity, we calculate thousands of random predictions together with the corresponding GALPROP results for the mentioned CR fluxes, and show the 1$\sigma$ uncertainty band of their relative differences in comparison with the AMS-02 measurement uncertainties in Fig. 2.

For both the proton and secondary antiproton fluxes, the differences between the prediction and GALPROP results are always smaller than the measurement uncertainties by one order of magnitude or more. The tertiary antiprotons that are produced by the interaction between the secondary antiprotons and ISM hold larger relative differences than the other two cases. Fortunately, these differences would not affect the observed antiproton flux, since the tertiary fluxes are much lower than the secondary fluxes.

Although only the training accuracies for the DR-2 + EPOS LHC model are shown in this subsection, we have ensured that all the models discussed in this work result in similar training accuracies.

III.3 Training for the DM contribution

The training procedure for the DM contribution is similar to that in Sec. III.2 except for two differences. Firstly, the DM contribution is irrelevant with the injection of CR particles, thus we only need to involve the propagation parameters in the training. Secondly, the shape of the DM injected spectrum depends on both the annihilation channel and DM mass $m_{\chi}$, thus we need to separately perform the training for each injection. In order to approach a more flexible predicting model, we just train a series of fluxes resulting from the injection of delta function in each energy bin, and then linearly combine them according to the DM injected spectrum. These injection models could only be used to predict the DM contribution, since their spatial distribution is set according to the NFW profile.

As a result, We show a random prediction example for the $b\bar{b}$ channel with a annihilating DM mass of $30\,\mathrm{GeV}$ in the DR-2 propagation model in Fig. 3. The DM contributions are much smaller than the measurements, thus we do not compare the corresponding relative difference with the measurement uncertainty. It can be seen that these predictions are even more accurate than the background predictions with relative differences at level of $0.01\%$. This is because that the predicting models for the DM contribution are simpler as they are constructed with only the propagation parameters.

IV Implications for DM annihilation

In this section we determine the confidence regions in the $m_{\chi}-\langle\sigma v\rangle$ plane favored by the AMS-02 antiproton measurement by scanning the parameter space. For the 2-dimensional parameter space of $(m_{\chi},\langle\sigma v\rangle)$, the boundary of the 1$\sigma$, 2$\sigma$ and 3$\sigma$ regions is given by varying the profile likelihoods as $\chi^{2}_{1\sigma,2\sigma,3\sigma}=\chi^{2}_{\mathrm{best}}+2.29$, 6.16, 11.6, respectively. We use the package NLOPT NLO to search for the best-fit and the required boundaries in the $m_{\chi}-\langle\sigma v\rangle$ plane.

In order to find the corresponding profile likelihood for each given point $(m_{\chi},\langle\sigma v\rangle)$, we optimize throughout the 2$\sigma$ confidence region of the astrophysical parameters constrained by the $\mathrm{B}/\mathrm{C}$, $\mathrm{Be}^{10}/\mathrm{Be}^{9}$ and proton data. Using the predicting models built in the above section, we could determine whether a set of propagation + injection + modulation parameters $\bm{\theta}$ falls inside the 2$\sigma$ confidence region, and reproduce the antiproton flux accurately within a time scale of $\mathcal{O}(0.1)$s. This optimization limited in the constrained astrophysical parameter region is also approached using the NLOPT with the augmented Lagrangian algorithm Birgin and Martínez (2008). The modulation potential for negative charged particles $\phi$ is allowed to vary from 0.5 to 1.5 times of $\phi_{+}$ as mentioned in Sec. II.2. Since the uncertainties associated with the heavy nuclear reaction cross sections, such as $\mathrm{C,O\rightarrow B}$, would affect the $\mathrm{B}/\mathrm{C}$ prediction and the determination of the propagation parameters, we also introduce a free scale factor $c_{\bar{p}}$ for the antiproton flux following Ref. Lin et al. (2016).

We investigate the DM implication in all the chosen hadronic models and propagation models. The results for the empirical and phenomenological hadronic models are shown in Fig. 4 and Fig. 5, respectively. We consider both the DR and DR-2 propagation models in these figures. The different colors in the figures represent the $1\sigma,2\sigma,3\sigma$ regions favored by the AMS-02 data. These regions that are not closed but extend to zero indicate upper limits at the corresponding significance instead of the interpretation of DM signatures.

We find that the differences arising from hadronic models are always larger than those from propagation models. It can be clearly seen from these figures that in most cases the significance of a DM signal is smaller than $2\sigma$. Among all the hadronic models, only the EPOS LHC model requires a DM contribution with $m_{\chi}\sim 1\,\mathrm{TeV}$ and $\langle\sigma v\rangle\sim 5\times 10^{25}\,\mathrm{cm}^{3}/\,\mathrm{s}$ at a confidence level over $99.7\%$, and result in an enclosed 3$\sigma$ confidence region. Unfortunately the whole 3$\sigma$ confidence region is excluded by the combined limit of the dwarf galaxy $\gamma$-ray emission set by Fermi and Magic Ahnen et al. (2016). However, the upper limit on the DM annihilation rate by the $\gamma$-ray emission from dwarf galaxies is model dependent, such as the case discussed in Ref. Xiang et al. (2017). If the DM annihilation is velocity dependent, it is still possible to contribute to the AMS-02 antiproton flux. The significance of the DM contribution for the other models is less than 2$\sigma$. For the Winkler’s model and the EPOS 1.99 model, the confidence regions expand throughout the lower panel and their upper boundaries become upper-limits on the DM annihilation rate.

A bump around $100\,\mathrm{GeV}$ is found in the upper boundaries for the EPOS 1.99 and Winkler’s model due to the underestimation of the secondary CR antiproton below $\sim 10\,\mathrm{GeV}$. This result is different from Refs. Cui et al. (2018b); Cuoco et al. (2019) in which the underestimation is quite significant and an enclosed contour instead of a bump is found. The reason is that those two works take the DR propagation model and adopt a charge independent solar modulation potential.

For comparison, we also show the confidence regions derived under the assumption $\phi_{-}=\phi_{+}$ for the Winkler’s model in Fig. 6. It can be seen that the confidence region for the DR-2 model stays unchanged while that for the DR model becomes an enclosed contour around $100\,\mathrm{GeV}$. The DR model tends to underestimate the secondary CR antiprotons at low energies, which is a well known fact and has been discussed very early Moskalenko et al. (2003); Trotta et al. (2011); Hooper et al. (2015). The deficiency of antiprotons at low energy is one of the main reasons to propose the DR-2 model Evoli et al. (2012). Our study shows that this deficiency can be corrected by a charge dependent solar modulation. A question is whether our allowing range for the negative charged potential ($0.5\phi_{+}\sim 1.5\phi_{+}$) is suitable for the charge-sign dependent modulation effect. Different simplifications on the solar modulation do affect the implication of the DR models. For instance, adopting a parametrized modulation potential to describe the charge-sign dependent effect, Cholis et. al. have found that the deficiency is still significant Cholis et al. (2019). A thorough solution to this problem is numerically solving the transport equation in the Solar system Maccione (2013); Potgieter (2014); Kappl (2016) instead of the FFA method. Such study is out of scope of this paper and is left for the future work.

In all, the significance of the DM contribution in the antiproton flux depends on the hadronic models. Only the EPOS LHC model would significantly expect the DM contribution at high energies. Note that the hadronic models are being developed with the collection of new collider data. For instance, a measurement of the antiproton production in $p-\mathrm{He}$ collisions at $\sqrt{s}=110\,\mathrm{GeV}$ has been recently presented by the LHCb collaboration Aaij et al. (2018), which shows a deviation from the predictions of many current hadronic models. Korsmeier et al. have developed the model of di Mauro et al. to fit this measurement Korsmeier et al. (2018). In the future, more relevant collider data will be collected and the hadronic models discussed above would be updated. The implication for the DM contribution in the CR antiproton flux will become more clear.

V Conclusions

In this work, we used the SVR method to rapidly predicted the CR antiproton fluxes with either the astrophysical or DM origination, and adopted them to investigate the DM parameter regions favored by the AMS-02 antiproton measurement. In the implementation, we have ensured that these predicted antiproton fluxes are accurate approximations of the GALPROP solutions in comparison with the current measurement uncertainties. The uncertainties of the astrophysical parameters have been well taken into account in the scanning analysis.

The DM implications in several propagation and hadronic models have been separately investigated. These models would occasionally underestimate the antiproton fluxes below several $\,\mathrm{GeV}$ or around $100\,\mathrm{GeV}$, and leave the room for the DM contribution with a DM mass around $100\,\mathrm{GeV}$ or $1\,\mathrm{TeV}$, respectively. If the charge-sign dependent solar modulation is taken into account, the favored region for a DM mass around $\sim 100\,\mathrm{GeV}$ could become not significant in comparison with the results in Refs. Cui et al. (2018b); Cuoco et al. (2019). However, this conclusion depends on the simplified charge-sign dependent solar modulation model used in our analysis. Further study based on the numerical solar modulation solution is needed to verify it. On the other hand, the implication for the $\,\mathrm{TeV}$ DM depends on the hadronic models. The DM contribution with a mass of $\,\mathrm{TeV}$ is only favored in the EPOS LHC model with a significance over 3$\sigma$. The future collider measurement will help to update these hadronic models and leave us a more certain conclusion on the DM implication in the CR antiprotons.

Acknowledgements.

This work is supported by the National Key R&D Program of China (No. 2016YFA0400200), the National Natural Science Foundation of China (Nos. U1738209 and 11851303).

Appendix A Support-vector machine

The SVM algorithm was originally invented by V. N. Vapnik and A. Y. Chervonenkis in 1963 as a binary linear classifier. Given a dataset of $n$ points in $l$ dimensional space $\{(\mathbf{x}_{i},y_{i})|i=1,2\dots n\}$ where the $y_{i}$s are either 1 or -1, the original SVM algorithm would tend to find two parallel boundary hyperplanes that could separate the two kinds of points with a maximum margin. These points are not promised to be linearly separable, thus a “soft-margin” rather than the real margin is adopted in the searching. The algorithm could then mathematically expressed as

[TABLE]

where $\lambda$ is a positive number determined manually. When the points are linearly separable, it is able to find a suitable $\bm{\omega}$ that keep all the $\max(0,1-y_{i}(\bm{\omega}\cdot\bm{x}_{i}-b))$ equal to 0, and $\bm{\omega}$ at the infimum in Eq. 6 naturally corresponds to the maximum margin. When they are not separable, the sum would always be positive, and $\lambda$ is used to trade-off between the classification error and the margin size. At the final solution, only a part of points that satisfy $1-y_{i}(\bm{\omega}\cdot\bm{x}_{i}-b)>0$ would directly contribute to the value of $\bm{\omega}$ and $b$. These points are called as ”support-vector”.

Based on the SVM algorithm, the SVR algorithm was then introduced in 1996 Drucker et al. (1996)

[TABLE]

where the value of $y_{i}$ is no longer dispersed, and the positive $\varepsilon$ is introduced to control the tolerance between $y_{i}$ and the predicted value $\bm{\omega}\cdot\bm{x}_{i}+b$. Smaller $\varepsilon$ would lead to more support-vectors and vise versa.

Eq. 7 is a linear regression; in order to deal with the non-linear problems, a kernel trick is introduced Boser et al. (1992). Any non-linear function $\hat{y}(\bm{x})$ could be transformed to a linear one by projecting the input parameter $\bm{x}$ to a space with higher dimension using a transform $\hat{\bm{\varphi}}(\bm{x})$, such that

[TABLE]

The kernel function is defined by the transform $\hat{\bm{\varphi}}(\bm{x})$ as

[TABLE]

In the implementation, we could directly adopt the kernel function to avoid the complicated calculation of $\hat{\bm{\varphi}}(\bm{x})$ when calculating the dot product in the transformed space $\mathcal{R}^{k}$. In this work, we adopt the widely used Gaussian radial basis function (RBF)

[TABLE]

where $\gamma$ is a positive free parameter. The RBF corresponds to a transformed space with infinite dimension Shashua (2009), thus it could be used to deal with almost all kinds of target functions.

Therefore, the SVR for non-linear problem is

[TABLE]

When the dataset size $n$ is large, directly searching for the infimum in Eq. 11 would require us to calculate the kernel function for many times. In order to avoid the time-consuming calculation for the kernel function, the problem is always transformed to its Lagrangian dual before being solved. We can express Eq. 11 as

[TABLE]

The corresponding Lagrangian is then

[TABLE]

where $\bm{W}=\hat{\bm{\varphi}}(\bm{\omega})$. The minimization that we search for is just the saddle point of $\mathcal{L}(\bm{W},b,\xi^{*},\xi,\alpha^{*},\alpha,\gamma)$ with respect to $\bm{W}$, $b$, $\xi^{*}$ and $\xi$, in which $\partial\mathcal{L}/\partial\bm{W}=\partial\mathcal{L}/\partial b=\partial\mathcal{L}/\partial\xi=\partial\mathcal{L}/\partial\xi^{*}=0$ Shashua (2009). Substituting these expressions back into Eq. 13, we could finally reach the dual problem

[TABLE]

The approximate function for $\hat{y}(\bm{x})$ is then $\hat{h}(\bm{x})=\sum_{i=1}^{n}(\alpha^{*}_{i}-\alpha_{i})k(\bm{x}_{i},\bm{x})+b$. This dual problem is a simple quadratic programming problem that could be solved fast with the sequential minimal optimization (SMO) Platt (1998).