Scrutinizing the evidence for dark matter in cosmic-ray antiprotons

TL;DR

This study critically examines the robustness of potential dark matter signals in cosmic-ray antiprotons by analyzing uncertainties in production cross sections and data correlations, finding that data correlations significantly impact the significance of the signal.

Contribution

The paper systematically evaluates the effects of cross-section uncertainties and data correlations on dark matter signal significance in cosmic-ray antiproton data, highlighting the importance of covariance information.

Findings

Cross-section uncertainties have a small effect on the dark matter signal significance.

Including data correlations can increase the signal significance from 3σ to above 5σ.

Proper covariance data is crucial for accurate interpretation of cosmic-ray measurements.

Abstract

Global fits of primary and secondary cosmic-ray (CR) fluxes measured by AMS-02 have great potential to study CR propagation models and search for exotic sources of antimatter such as annihilating dark matter (DM). Previous studies of AMS-02 antiprotons revealed a possible hint for a DM signal which, however, could be affected by systematic uncertainties. To test the robustness of such a DM signal, in this work we systematically study two important sources of uncertainties: the antiproton production cross sections needed to calculate the source spectra of secondary antiprotons and the potential correlations in the experimental data, so far not provided by the AMS-02 Collaboration. To investigate the impact of cross-section uncertainties we perform global fits of CR spectra including a covariance matrix determined from nuclear cross-section measurements. As an alternative approach, we perform a joint fit to both the CR and cross-section data. The two methods agree and show that cross-section uncertainties have a small effect on the CR fits and on the significance of a potential DM signal, which we find to be at the level of $3\sigma$. Correlations in the data can have a much larger impact. To illustrate this effect, we determine possible benchmark models for the correlations in a data-driven method. The inclusion of correlations strongly improves the constraints on the propagation model and, furthermore, enhances the significance of the DM signal up to above $5\sigma$. Our analysis demonstrates the importance of providing the covariance of the experimental data, which is needed to fully exploit their potential.

Figures

Tables

Table 1. Table 1: Fit quality for the best-fit parameter points without (second column) and with (third column) DM for various choices of the rigidity cut. The last column shows the absolute Δ ​ χ 2 Δ superscript 𝜒 2 \Delta\chi^{2} between the respective fits with and without DM.

	${\chi }^2/\text{ndf}$
rigidity cut [GV]	excl. DM	incl. DM	$\Delta {\chi }^2$	DM significance
5	$35.6/145=0.245$	$22.9/143=0.160$	12.7	$3.1\sigma$
3	$52.7/160=0.329$	$34.2/158=0.216$	18.5	$3.9\sigma$
2	$68.2/172=0.396$	$57.1/170=0.336$	11.1	$2.9\sigma$
1	$105.4/182=0.579$	$105.6/180=0.586$	-0.2	–

Table 2. Table 2: Covariance matrices for the different benchmark scenarios. In the case of proton and helium, we maximize the log-likelihood in Eq. ( 22 ) at fixed ℓ corr subscript ℓ corr \ell_{\mathrm{corr}} to determine f 𝑓 f and α 𝛼 \alpha , while, to be conservative, we fix f p ¯ subscript 𝑓 ¯ 𝑝 f_{\bar{p}} to 0.3 and α p ¯ subscript 𝛼 ¯ 𝑝 \alpha_{\bar{p}} to 1.

		${\ell }_{\text{corr}}=\mathrm{𝟎}$					${\ell }_{\text{corr}}=\mathrm{𝟓}$					${\ell }_{\text{corr}}=\mathrm{𝟏𝟎}$
dataset		$p$	He	$\overline{p}/p$			$p$	He	$\overline{p}/p$			$p$	He	$\overline{p}/p$
$f$		0.062	0.080	0.30			0.079	0.103	0.30			0.082	0.101	0.30
$\alpha$		0.63	0.81	1.00			0.20	0.21	1.00			0.20	0.20	1.00

Table 3. Table 3: The best-fit parameters of the fits with and without DM. The double-column contains the results of our default setup. We show for comparison a fit that includes CREAM data and AMS-02 above 300 GeV.

parameter	default setup		extended setup
XS parametrization	param. MW		param. MW
DM	incl.	excl.	incl.	excl.
${\gamma }_1$	${1.71}_{-0.25}^{+0.02}$	${1.72}_{-0.12}^{+0.04}$	${1.64}_{-0.12}^{+0.05}$	${1.72}_{-0.11}^{+0.02}$
${\gamma }_{1,p}$	${1.78}_{-0.19}^{+0.003}$	${1.75}_{-0.10}^{+0.03}$	${1.73}_{-0.06}^{+0.04}$	${1.73}_{-0.07}^{+0.05}$
${\gamma }_2$	${2.41}_{-0.002}^{+0.03}$	${2.38}_{-0.02}^{+0.01}$	${2.44}_{-0.02}^{+0.01}$	${2.38}_{-0.01}^{+0.01}$
${\gamma }_{2,p}$	${2.45}_{-0.002}^{+0.03}$	${2.42}_{-0.02}^{+0.01}$	${2.48}_{-0.02}^{+0.01}$	${2.41}_{-0.01}^{+0.01}$
$R_0,\left[\mathrm{MV}\right]$	${6950}_{-1640}^{+330}$	${7380}_{-1450}^{+910}$	${6519}_{-824}^{+1045}$	${7695}_{-1375}^{+563}$
$s_0$	${0.38}_{-0.04}^{+0.06}$	${0.34}_{-0.04}^{+0.05}$	${0.38}_{-0.01}^{+0.06}$	${0.37}_{-0.03}^{+0.04}$
$D_0,\left[{10}^{28}{\mathrm{cm}}^2/s\right]$	${5.43}_{-3.17}^{+0.45}$	${2.90}_{-1.21}^{+1.33}$	${3.46}_{-1.19}^{+2.46}$	${1.97}_{-3.81}^{+1.27}$
$\delta$	${0.38}_{-0.03}^{+0.01}$	${0.42}_{-0.01}^{+0.02}$	${0.35}_{-0.01}^{+0.02}$	${0.42}_{-0.02}^{+0.01}$
$v_A,\left[\mathrm{km}/h\right]$	${18.0}_{-1.4}^{+2.1}$	${16.2}_{-2.5}^{+1.0}$	${18.59}_{-3.25}^{+0.00}$	${15.81}_{-1.99}^{+0.87}$
$v_{0,c},\left[\mathrm{km}/h\right]$	${0.08}_{-0.08}^{+9.09}$	${0.52}_{-0.51}^{+2.32}$	${0.35}_{-0.14}^{+4.94}$	${0.79}_{-0.77}^{+2.19}$
$z_h,\left[\mathrm{kpc}\right]$	${6.45}_{-4.26}^{+0.30}$	${3.58}_{-1.52}^{+2.36}$	${3.36}_{-1.13}^{+3.47}$	${2.47}_{-0.43}^{+1.58}$
$\log \left(m_{\mathrm{DM}}/\left[\mathrm{GeV}\right]\right)$	${1.89}_{-0.08}^{+0.03}$		${1.88}_{-0.00}^{+0.05}$
$\log \left(⟨\sigma v⟩/\left[s/{\mathrm{cm}}^3\right]\right)$	$-{26.16}_{-0.04}^{+0.78}$		$-{25.56}_{-0.47}^{+0.20}$
${\delta }_2$			${0.23}_{-0.00}^{+0.01}$	${0.30}_{-0.02}^{+0.02}$
$R_1,\left[\mathrm{GV}\right]$			${344}_{-20}^{+26}$	${338}_{-40}^{+49}$
${\phi }_{\text{SM,AMS-02},p,\mathrm{He}},\left[\mathrm{MV}\right]$	${616}_{-72}^{+71}$	${625}_{-85}^{+55}$	${566}_{-71}^{+19}$	${567}_{-66}^{+19}$
${\phi }_{\text{SM,AMS-02},\overline{p}},\left[\mathrm{MV}\right]$	${604}_{-114}^{+112}$	${561}_{-112}^{+135}$	${577}_{-82}^{+43}$	${561}_{-106}^{+35}$
${\chi }_{\text{AMS-02},p}^2$	3.2	2.6	4.9	7.0
${\chi }_{\text{AMS-02,He}}^2$	4.0	4.8	10.4	12.0
${\chi }_{\text{AMS-02},\overline{p}}^2$	11.1	22.1	12.2	20.0
${\chi }_{\text{Voager},p}^2$	3.2	3.8	2.7	3.6
${\chi }_{\text{Voager,He}}^2$	1.3	1.9	3.4	4.1
${\chi }_{\text{CREAM},p}^2$			1.0	1.0
${\chi }_{\text{CREAM,He}}^2$			1.0	2.9
${\chi }_{{\phi }_{\text{SM}}}^2$	0.0	0.4	0.0	0.0
${\chi }_{\text{CR}}^2$	22.9	35.6	35.6	50.7
${\chi }^2$/dof	22.9/143	35.6/145	35.6/177	50.7/179
$\Delta {\chi }^2$	12.7		15.1
DM significance	$3.1\sigma$		$3.5\sigma$

Table 4. Table 4: The best-fit parameters of various fits to test the impact of cross-section uncertainties. For details refer to Section  IV .

parameter	default setup, different XS		XS uncertainty by covariance matrix				joint fit (CR$+$XS)
XS parametrization	param. MD		param. MW		param. MD		param. MW
DM	incl.	excl.	incl.	excl.	incl.	excl.	incl.	excl.
${\gamma }_1$	${1.55}_{-0.05}^{+0.11}$	${1.65}_{-0.09}^{+0.09}$	${1.65}_{-0.21}^{+0.07}$	${1.68}_{-0.09}^{+0.07}$	${1.52}_{-0.06}^{+0.08}$	${1.64}_{-0.06}^{+0.09}$	${1.57}_{-0.05}^{+0.06}$	${1.67}_{-0.09}^{+0.09}$
${\gamma }_{1,p}$	${1.67}_{-0.02}^{+0.08}$	${1.69}_{-0.07}^{+0.09}$	${1.76}_{-0.17}^{+0.03}$	${1.73}_{-0.08}^{+0.05}$	${1.64}_{-0.05}^{+0.08}$	${1.69}_{-0.05}^{+0.07}$	${1.68}_{-0.05}^{+0.06}$	${1.70}_{-0.05}^{+0.07}$
${\gamma }_2$	${2.43}_{-0.01}^{+0.01}$	${2.39}_{-0.02}^{+0.01}$	${2.43}_{-0.02}^{+0.02}$	${2.37}_{-0.01}^{+0.02}$	${2.43}_{-0.01}^{+0.02}$	${2.28}_{-0.01}^{+0.01}$	${2.44}_{-0.03}^{+0.01}$	${2.38}_{-0.03}^{+0.01}$
${\gamma }_{2,p}$	${2.48}_{-0.01}^{+0.01}$	${2.42}_{-0.02}^{+0.01}$	${2.47}_{-0.02}^{+0.02}$	${2.40}_{-0.01}^{+0.02}$	${2.47}_{-0.01}^{+0.02}$	${2.41}_{-0.01}^{+0.01}$	${2.48}_{-0.03}^{+0.01}$	${2.42}_{-0.03}^{+0.01}$
$R_0,\left[\mathrm{MV}\right]$	${5860}_{-300}^{+870}$	${6700}_{-620}^{+1490}$	${6910}_{-1710}^{+470}$	${6880}_{-940}^{+1200}$	${5840}_{-430}^{+510}$	${6830}_{-130}^{+730}$	${5950}_{-370}^{+640}$	${6910}_{-890}^{+970}$
$s_0$	${0.43}_{-0.08}^{+0.01}$	${0.41}_{-0.06}^{+0.02}$	${0.38}_{-0.02}^{+0.10}$	${0.36}_{-0.06}^{+0.04}$	${0.41}_{-0.03}^{+0.04}$	${0.38}_{-0.04}^{+0.02}$	${0.43}_{-0.08}^{+0.01}$	${0.37}_{-0.05}^{+0.04}$
$D_0,\left[{10}^{28}{\mathrm{cm}}^2/s\right]$	${2.70}_{-0.25}^{+2.72}$	${2.10}_{-0.50}^{+0.89}$	${5.64}_{-3.60}^{+0.50}$	${2.94}_{-1.06}^{+0.68}$	${2.27}_{-0.13}^{+0.90}$	${1.86}_{-0.19}^{+0.61}$	${2.47}_{-0.01}^{+4.21}$	${1.72}_{-0.26}^{+1.81}$
$\delta$	${0.35}_{-0.01}^{+0.02}$	${0.41}_{-0.02}^{+0.02}$	${0.35}_{-0.02}^{+0.03}$	${0.43}_{-0.03}^{+0.01}$	${0.35}_{-0.02}^{+0.02}$	${0.42}_{-0.01}^{+0.02}$	${0.35}_{-0.02}^{+0.04}$	${0.42}_{-0.0004}^{+0.02}$
$v_A,\left[\mathrm{km}/h\right]$	${16.8}_{-1.8}^{+2.7}$	${14.6}_{-1.7}^{+1.9}$	${19.5}_{-4.9}^{+0.5}$	${15.0}_{-0.6}^{+2.1}$	${17.2}_{-2.2}^{+2.7}$	${14.3}_{-1.8}^{+1.4}$	${17.9}_{-2.2}^{+0.3}$	${15.2}_{-1.38}^{+0.66}$
$v_{0,c},\left[\mathrm{km}/h\right]$	${4.80}_{-4.41}^{+2.44}$	${1.97}_{-1.95}^{+4.11}$	${2.63}_{-2.44}^{+6.21}$	${0.66}_{-0.61}^{+1.65}$	${4.22}_{-3.66}^{+3.67}$	${2.78}_{-2.40}^{+1.01}$	${1.17}_{-0.15}^{+5.39}$	${0.27}_{-0.19}^{+2.02}$
$z_h,\left[\mathrm{kpc}\right]$	${2.62}_{-0.23}^{+2.49}$	${2.54}_{-0.53}^{+1.24}$	${5.76}_{-3.70}^{+0.50}$	${3.72}_{-1.48}^{+0.82}$	${2.20}_{-0.14}^{+0.81}$	${2.28}_{-0.22}^{+0.86}$	${2.46}_{-0.11}^{+4.21}$	${2.15}_{-0.15}^{+2.61}$
$\log \left(m_{\mathrm{DM}}/\left[\mathrm{GeV}\right]\right)$	${1.85}_{-0.05}^{+0.02}$		${1.88}_{-0.06}^{+0.07}$		${1.82}_{-0.03}^{+0.04}$		${1.89}_{-0.03}^{+0.04}$
$\log \left(⟨\sigma v⟩/\left[s/{\mathrm{cm}}^3\right]\right)$	$-{25.50}_{-0.36}^{+0.01}$		$-{25.88}_{-0.09}^{+0.59}$		$-{25.32}_{-0.41}^{+0.01}$		$-{25.47}_{-0.42}^{+0.10}$
${\phi }_{\text{SM,AMS-02},p,\mathrm{He}},\left[\mathrm{MV}\right]$	${554}_{-13}^{+101}$	${544}_{-58}^{+44}$	${582}_{-64}^{+90}$	${582}_{-25}^{+54}$	${582}_{-46}^{+38}$	${556}_{-28}^{+7}$	${584}_{-80}^{+39}$	${598}_{-61}^{+22}$
${\phi }_{\text{SM,AMS-02},\overline{p}},\left[\mathrm{MV}\right]$	${549}_{-45}^{+31}$	${522}_{-115}^{+62}$	${592}_{-90}^{+120}$	${605}_{-114}^{+56}$	${600}_{-56}^{+54}$	${553}_{-55}^{+17}$	${498}_{-81}^{+166}$	${537}_{-61}^{+152}$
$C_1,\left[{10}^{-3}{\left(\mathrm{GeV}\right)}^{-2}\right]$							${50.1}_{-1.0}^{+0.5}$	${51.0}_{-1.7}^{+0.5}$
$C_5,\left[{10}^{-3}\right]$							${0.32}_{-0.02}^{+0.24}$	${0.27}_{-0.17}^{+0.08}$
$C_6$							${3.82}_{-0.40}^{+0.27}$	${3.73}_{-0.10}^{+0.81}$
${\chi }_{\text{AMS-02},p}^2$	2.4	3.8	2.6	2.6	2.2	2.6	2.8	2.6
${\chi }_{\text{AMS-02,He}}^2$	4.4	4.6	4.1	5.1	4.1	4.3	3.9	4.7
${\chi }_{\text{AMS-02},\overline{p}}^2$	11.0	21.1	11.4	21.3	11.3	22.4	11.7	19.9
${\chi }_{\text{Voager},p}^2$	3.5	3.7	3.2	3.7	3.1	3.8	2.8	3.5
${\chi }_{\text{Voager,He}}^2$	0.9	1.0	1.7	1.1	1.1	1.0	1.1	1.9
${\chi }_{{\phi }_{\text{SM}}}^2$	0.0	0.1	0.0	0.1	0.0	0.0	0.7	0.4
${\chi }_{\text{CR}}^2$	22.2	34.2	23.0	33.9	22.0	34.2	23.2	33.0
${\chi }_{\text{XS}}^2$							791.4	792.2
${\chi }^2$/dof	22.2/143	34.2/145	23.0/143	33.9/145	22.0/143	34.2/145	814.5/799	825.2/801
$\Delta {\chi }^2$	12.0		10.9		12.2		10.7
DM significance	$3.0\sigma$		$2.9\sigma$		$3.1\sigma$		$2.8\sigma$

Table 5. Table 5: The best-fit parameters of various fits to test the impact of correlated AMS-02 uncertainties. For details refer to Section  V .

parameter	1%${\sigma }_{\text{sys}}$		${\ell }_{\text{corr}}=\mathrm{𝟎}$		${\ell }_{\text{corr}}=\mathrm{𝟓}$		${\ell }_{\text{corr}}=\mathrm{𝟏𝟎}$
XS parametrization	param. MW		param. MW		param. MW		param. MW
DM	incl.	excl.	incl.	excl.	incl.	excl.	incl.	excl.
${\gamma }_1$	${1.62}_{-0.10}^{+0.08}$	${1.77}_{-0.08}^{+0.02}$	${1.65}_{-0.03}^{+0.02}$	${1.70}_{-0.05}^{+0.03}$	${1.61}_{-0.12}^{+0.10}$	${1.72}_{-0.06}^{+0.02}$	${1.65}_{-0.11}^{+0.07}$	${1.67}_{-0.06}^{+0.03}$
${\gamma }_{1,p}$	${1.70}_{-0.07}^{+0.08}$	${1.74}_{-0.07}^{+0.02}$	${1.67}_{-0.04}^{+0.01}$	${1.66}_{-0.04}^{+0.02}$	${1.64}_{-0.13}^{+0.09}$	${1.69}_{-0.06}^{+0.02}$	${1.66}_{-0.13}^{+0.04}$	${1.65}_{-0.06}^{+0.02}$
${\gamma }_2$	${2.43}_{-0.01}^{+0.01}$	${2.37}_{-0.002}^{+0.02}$	${2.42}_{-0.01}^{+0.00}$	${2.39}_{-0.00}^{+0.01}$	${2.43}_{-0.01}^{+0.00}$	${2.40}_{-0.00}^{+0.00}$	${2.42}_{-0.01}^{+0.01}$	${2.39}_{-0.00}^{+0.01}$
${\gamma }_{2,p}$	${2.47}_{-0.01}^{+0.01}$	${2.40}_{-0.001}^{+0.02}$	${2.45}_{-0.01}^{+0.00}$	${2.41}_{-0.00}^{+0.01}$	${2.46}_{-0.02}^{+0.00}$	${2.42}_{-0.00}^{+0.00}$	${2.45}_{-0.01}^{+0.01}$	${2.42}_{-0.00}^{+0.01}$
$R_0,\left[\mathrm{MV}\right]$	${6600}_{-1090}^{+720}$	${9180}_{-1270}^{+50}$	${6869}_{-173}^{+200}$	${7844}_{-512}^{+575}$	${6519}_{-793}^{+1300}$	${7983}_{-707}^{+340}$	${6743}_{-743}^{+1103}$	${7290}_{-582}^{+262}$
$s_0$	${0.42}_{-0.05}^{+0.03}$	${0.32}_{-0.01}^{+0.06}$	${0.45}_{-0.01}^{+0.01}$	${0.40}_{-0.01}^{+0.02}$	${0.47}_{-0.04}^{+0.01}$	${0.41}_{-0.01}^{+0.02}$	${0.46}_{-0.03}^{+0.02}$	${0.41}_{-0.01}^{+0.02}$
$D_0,\left[{10}^{28}{\mathrm{cm}}^2/s\right]$	${2.80}_{-0.85}^{+2.35}$	${2.19}_{-0.58}^{+0.36}$	${1.99}_{-0.02}^{+0.05}$	${1.79}_{-0.013}^{+0.34}$	${1.94}_{-0.20}^{+0.22}$	${1.90}_{-0.03}^{+0.05}$	${1.84}_{-0.02}^{+0.22}$	${1.66}_{-0.00}^{+0.25}$
$\delta$	${0.36}_{-0.01}^{+0.02}$	${0.43}_{-0.02}^{+0.002}$	${0.37}_{-0.01}^{+0.01}$	${0.42}_{-0.01}^{+0.00}$	${0.36}_{-0.00}^{+0.02}$	${0.40}_{-0.00}^{+0.00}$	${0.37}_{-0.01}^{+0.01}$	${0.41}_{-0.01}^{+0.00}$
$v_A,\left[\mathrm{km}/h\right]$	${21.1}_{-3.4}^{+0.4}$	${19.6}_{-1.4}^{+1.2}$	${21.87}_{-0.17}^{+0.17}$	${20.93}_{-0.45}^{+1.50}$	${22.23}_{-1.77}^{+1.89}$	${22.15}_{-0.93}^{+0.15}$	${21.37}_{-0.79}^{+1.51}$	${20.25}_{-0.17}^{+1.41}$
$v_{0,c},\left[\mathrm{km}/h\right]$	${0.52}_{-0.27}^{+4.62}$	${0.06}_{-1.36}^{+0.79}$	${0.53}_{-0.27}^{+2.75}$	${0.37}_{-0.26}^{+1.51}$	${3.69}_{-3.66}^{+10.95}$	${0.24}_{-0.21}^{+2.65}$	${2.01}_{-1.92}^{+9.99}$	${0.25}_{-0.00}^{+3.96}$
$z_h,\left[\mathrm{kpc}\right]$	${2.76}_{-0.75}^{+2.68}$	${2.72}_{-0.70}^{+0.44}$	${2.03}_{-0.00}^{+0.08}$	${2.05}_{-0.02}^{+0.41}$	${2.00}_{-0.00}^{+0.29}$	${2.11}_{-0.01}^{+0.03}$	${2.01}_{-0.00}^{+0.12}$	${2.02}_{-0.02}^{+0.13}$
$\log \left(m_{\mathrm{DM}}/\left[\mathrm{GeV}\right]\right)$	${1.82}_{-0.02}^{+0.05}$		${1.70}_{-0.01}^{+0.04}$		${1.74}_{-0.06}^{+0.01}$		${1.71}_{-0.06}^{+0.05}$
$\log \left(⟨\sigma v⟩/\left[s/{\mathrm{cm}}^3\right]\right)$	$-{25.57}_{-0.37}^{+0.26}$		$-{25.49}_{-0.12}^{+0.07}$		$-{25.41}_{-0.10}^{+0.10}$		$-{25.38}_{-0.14}^{+0.05}$
${\phi }_{\text{SM,AMS-02},p,\mathrm{He}},\left[\mathrm{MV}\right]$	${674}_{-138}^{+9}$	${738}_{-39}^{+58}$	${819}_{-0}^{+14}$	${857}_{-25}^{+16}$	${839}_{-43}^{+20}$	${844}_{-26}^{+14}$	${837}_{-12}^{+18}$	${832}_{-2}^{+26}$
${\phi }_{\text{SM,AMS-02},\overline{p}},\left[\mathrm{MV}\right]$	${601}_{-138}^{+146}$	${613}_{-94}^{+69}$	${741}_{-81}^{+71}$	${568}_{-55}^{+30}$	${743}_{-43}^{+190}$	${534}_{-26}^{+57}$	${791}_{-99}^{+92}$	${556}_{-34}^{+42}$
${\chi }_{\text{AMS-02},p}^2$	3.3	7.6	83.5	81.0	74.4	75.9	79.4	84.3
${\chi }_{\text{AMS-02,He}}^2$	7.6	9.8	97.4	98.9	96.5	94.5	96.7	89.7
${\chi }_{\text{AMS-02},\overline{p}}^2$	30.3	44.0	52.2	69.2	42.3	50.1	47.4	53.9
${\chi }_{\text{Voager},p}^2$	3.4	7.5	8.8	14.7	7.7	9.5	6.8	13.4
${\chi }_{\text{Voager,He}}^2$	2.1	6.8	7.3	11.4	9.3	10.1	10.4	10.0
${\chi }_{{\phi }_{\text{SM}}}^2$	0.6	1.7	0.9	8.8	2.4	10.0	0.6	8.0
${\chi }^2$/dof	47.4/145	77.4/143	250.0/143	284.1/145	232.6/143	250.2/145	241.3/143	259.3/145
$\Delta {\chi }^2$	30.0		34.1		17.6		18.0
DM significance	$5.1\sigma$		$5.5\sigma$		$3.8\sigma$		$3.8\sigma$

Equations

$$\frac{\partial {\psi }_i(\mathbf{x},p,t)}{\partial t}=q_i(\mathbf{x},p)$$

$$q_i(\mathbf{x},p)=q_i(r,z,R)=q_{0,i}{q}_{r,z}(r,z)q_{i,R}(R),$$

$$q_R(R)$$

$$q_{r,z}(r,z)={\left(\frac{r}{r_s}\right)}^{\alpha }\exp \left(-\beta \frac{r-r_s}{r_s}\right)\exp \left(-\frac{|z|}{z_0}\right),$$

$$D_{xx}$$

$$D_{pp}=\frac{4{\left(p{v}_{\mathrm{A}}\right)}^2}{3(2-\delta )(2+\delta )(4-\delta )\delta D_{xx}}.$$

$$q_{\bar{p}}^{(\mathrm{DM})}(\mathbf{x},T)=\frac{1}{2}{\left(\frac{\rho (\mathbf{x})}{m_{\mathrm{DM}}}\right)}^2\sum _f{⟨\sigma v⟩}_f\frac{\mathrm{d}{N}_{\bar{p}}^f}{\mathrm{d}T},$$

$${\varphi }_{\oplus ,i}(E_{\oplus ,i})$$

$$E_{\oplus ,i}$$

$$-2\log ({\mathcal{L}}_{\mathrm{CR}})={\chi }_{\mathrm{CR}}^2=\sum _{e,s}\sum _{i,j}\left({\varphi }_{s,i}^{(e)}-{\varphi }_{s,e}^{(\text{m})}(R_i)\right){\left({\left({\mathcal{V}}^{\left(e,s\right)}\right)}^{-1}\right)}_{ij}\left({\varphi }_{s,j}^{(e)}-{\varphi }_{s,e}^{(\text{m})}(R_j)\right).$$

$$q_{ij}(\mathbf{x},T_{\bar{p}})$$

$${\mathcal{V}}_{\text{XS},ij}^{(q_{\bar{p},\text{rel}})}=\frac{1}{N-1}\sum _{k=1}^N\left(\frac{q_{\bar{p}}^{(k)}(R_i)}{q_{\bar{p}}^{(\text{bestfit)}}(R_i)}-1\right)\left(\frac{q_{\bar{p}}^{(k)}(R_j)}{q_{\bar{p}}^{(\text{bestfit)}}(R_j)}-1\right).$$

$${\mathcal{V}}_{\text{XS},ij}^{\left({\varphi }_{\bar{p}/p}^{\text{AMS-02}}\right)}={\mathcal{V}}_{\text{XS},ij}^{(q_{\bar{p},\text{rel}})}{\varphi }_{\bar{p}/p,i}^{(\text{AMS-02})}{\varphi }_{\bar{p}/p,j}^{(\text{AMS-02})}.$$

$${\mathcal{V}}_{ij}^{\left({\varphi }_{\bar{p}/p}^{\text{AMS-02}}\right)}={\mathcal{V}}_{\text{XS},ij}^{\left({\varphi }_{\bar{p}/p}^{\text{AMS-02}}\right)}+{\delta }_{ij}{\left[{\sigma }_{\bar{p}/p,i}^{(\text{AMS-02})}\right]}^2.$$

$$E_{\bar{p}}\frac{{\mathrm{d}}^3\sigma }{\mathrm{d}{p}_{\bar{p}}^3}(\sqrt{s},x_{\mathrm{R}},p_{\mathrm{T}})$$

$$R_{\mathrm{XS}}$$

$$X$$

$$\log ({\mathcal{L}}_{\text{joint}})=\log ({\mathcal{L}}_{\mathrm{CR},\mathrm{SM}})+\log ({\mathcal{L}}_{\mathrm{XS}}).$$

$$-2\log ({\mathcal{L}}_{\mathrm{XS}})=\sum _e\sum _i{\left(\frac{{\omega }_e{\sigma }_{\text{inv},i}^{(e)}-{\sigma }_{\text{inv}}^{(m)}({\sqrt{s}}_i,{x_{\mathrm{R}}}_i,{p_{\mathrm{T}}}_i)}{{\omega }_e{\sigma }_{{\sigma }_{\text{inv},i}^{(e)}}}\right)}^2+\sum _e{\left(\frac{1-{\omega }_e}{{\sigma }_{{\omega }_e}}\right)}^2.$$

$$\mathcal{V}={\mathcal{V}}_{\text{stat}}+{\mathcal{V}}_{\text{short}}+{\mathcal{V}}_{\text{long}}.$$

$${\mathcal{V}}_{\text{short},ij}=\exp \left(-\frac{|i-j{|}^{\alpha }}{{{\ell }_{\text{corr}}}^{\alpha }}\right)f^2{\sigma }_{\text{sys},i}{\sigma }_{\text{sys},j},$$

$$-\log (\mathcal{L})=\frac{1}{2}\log \left(\det \left[{\mathcal{V}}_{\text{stat}}+{\mathcal{V}}_{\text{short}}({\ell }_{\text{corr}},f,\alpha )\right]\right)+\frac{1}{2}\mathbf{x}\cdot {\left[{\mathcal{V}}_{\text{stat}}+{\mathcal{V}}_{\text{short}}({\ell }_{\text{corr}},f,\alpha )\right]}^{-1}\cdot \mathbf{x}+\text{const},$$

Full text

††preprint: LAPTH-052/18, TTK-19-09, CP3-19-08

I Introduction

Antimatter in cosmic rays (CRs), and in particular antiprotons, have been extensively investigated as a powerful means to search for exotic CR sources, such as dark matter (DM) annihilation in the Galaxy Bergstrom et al. (1999); Donato et al. (2004); Bringmann and Salati (2007); Donato et al. (2009); Fornengo et al. (2014); Hooper et al. (2015); Pettorino et al. (2014); Boudaud et al. (2015); Cembranos et al. (2015); Cirelli et al. (2014); Bringmann et al. (2014); Giesen et al. (2015); Evoli et al. (2015, 2012); Reinert and Winkler (2018); Cui et al. (2018). The recent very accurate measurement of the CR antiproton flux by the AMS-02 experiment Aguilar et al. (2016) has significantly increased the sensitivity to a DM signal. A DM contribution as low as about $10\%$ of the antiproton flux can now in principle be detected, provided that the theoretical and experimental systematic uncertainties are under control at that level. Indeed, strong limits on heavy DM have been derived from global CR fits Cuoco et al. (2018). At the same time, the data have also revealed a tentative signal of DM, corresponding to a DM mass of around 40–130 GeV and a thermal annihilation cross section, $\left\langle\sigma v\right\rangle\sim 3\times 10^{-26}$ cm3/s Cuoco et al. (2017a); Cui et al. (2017); Cuoco et al. (2017b). This signal, if confirmed, is compatible with a DM interpretation of the Galactic center $\gamma$-ray excess (GCE) for a variety of annihilation channels. It is also expected to provide a detectable signal in antideuterons Korsmeier et al. (2018a), and it is compatible with a variety of different beyond-the-standard-model scenarios Clark et al. (2018); Cirelli et al. (2017); Cuoco et al. (2017b); Cline et al. (2017); Jia (2017); Arcadi et al. (2017); Li (2018); Escudero et al. (2017). However, given the small experimental errors, several important sources of systematic uncertainties, which before could be neglected, now become increasingly important and need to be further investigated.

One such uncertainty concerns the predictions for the antiproton cross sections, needed to model antiproton production through scattering of CR protons and helium with the interstellar medium (ISM) in the Galactic disk. Recent progress in the determination of antiproton cross sections from nuclear experimental data Winkler (2017) has been found to have a significant impact on the DM interpretation of cosmic rays Reinert and Winkler (2018). A second important source of uncertainty are possible correlations in the AMS-02 data, which are dominated by systematic uncertainties in most parts of the determined energy range.111The proton and helium fluxes are dominated by systematic uncertainties in the whole energy range from 1 GeV to 3 TeV, while in the antiproton-to-proton ratio systematic uncertainties are dominant from 1.8 to 50 GeV only. The AMS-02 Collaboration has only released absolute systematic uncertainties, without providing information about their correlations. However, most of the systematic uncertainties are expected to exhibit sizable correlations.

In this paper we systematically study the above-mentioned sources of uncertainties. We examine their impact on the significance of the potential DM excess in the CR antiproton data and hence scrutinize the robustness of this finding. Furthermore, we shed light on the impact of systematic uncertainties on the parameter determination of the CR propagation model. To this end, we first reproduce the finding of Cuoco et al. (2017a) in an updated setup, using the most recent cross-section parametrization from Korsmeier et al. (2018b) and an improved treatment of solar modulation. This fit is considered as the default setup providing the reference for the following investigations.

We study the impact of uncertainties in the antiproton production cross section following two approaches. The first approach is similar to the one taken in Reinert and Winkler (2018). We incorporate the cross-section uncertainties in the CR fit by including a covariance matrix extracted from a separate fit to nuclear measurements. In the second approach we perform a joint fit of the CR propagation and antiproton cross-section parameters to the AMS-02 and the nuclear data. Such a joint fit provides important information about possible correlation between propagation and cross-section parameters.

We investigate the potential impact of correlations in the experimental data by assuming that the systematic uncertainties consist of three components: a part which is uncorrelated, a part correlated over a certain number of neighboring rigidity bins, and a part which is fully correlated. We determine these properties in a data-driven method, which allows us to constrain the viable range of the various components and to define four corresponding benchmark models. We perform global fits for these four benchmark models, each with and without a primary source of antiprotons from DM.

The paper is structured as follows. In Sec. II we review CR propagation and highlight some new features of our setup which extend the standard treatment. Furthermore, we discuss solar modulation and justify our choice to omit low-rigidity data from the global fit. In Sec. III we detail the numerical implementation and describe our default setup for propagation and the corresponding results, considering the case with and without an additional source of antiprotons from DM annihilation. The uncertainties from antiproton cross sections are discussed in Sec. IV, following the two approaches mentioned above. Finally, in Sec. V we discuss the potential impact of correlated errors in the AMS-02 data, following again two different methods. We conclude in Sec. VI.

II Cosmic-ray propagation and solar modulation

The propagation of CRs through the Galaxy can be described by the well-known diffusion equation for the phase-space densities of CRs Strong et al. (2007):

[TABLE]

The equation can either be solved analytically, utilizing various simplifying assumption Putze et al. (2010); Maurin (2018), or fully numerically, as implemented in codes like Galprop Strong et al. (2000); Strong (2015), Dragon Evoli et al. (2008, 2018) and Picard Kissmann (2014).

For a given primary CR species $i$, the equation has a source term $q_{i}(\bm{x},p)$ which is assumed to factorize into a space- and rigidity-dependent part,

[TABLE]

where $r$ and $z$ are cylindrical coordinates with respect to the Galactic center and $R$ denotes the rigidity. The rigidity dependence is taken to be a smoothly broken power law with break position $R_{0}$, and spectral indices $\gamma_{1,i}$ and $\gamma_{2,i}$ above and below the break, respectively, while the smoothing is controlled by a parameter $s$:

[TABLE]

The spatial dependence of the source term is parametrized as

[TABLE]

with parameters $\alpha=0.5$, $\beta=2.2$, $r_{s}=8.5$ kpc, and $z_{0}=0.2$ kpc.222These are the default values in Galprop v56, which slightly differ from the values obtained from supernova remnants Case and Bhattacharya (1998); Green (2015). However, as pointed out in Korsmeier and Cuoco (2016) the spatial dependence of the source term distribution has a negligible effect on the local CR fluxes.

Several processes contribute to CR propagation, in particular diffusion, reacceleration, convection and energy losses. Spatial diffusion is assumed to be isotropic and homogeneous and is described by a rigidity-dependent diffusion coefficient $D_{xx}$. We use a broken power law in rigidity in light of the recent measurements by AMS-02 of the secondaries boron, lithium and beryllium Aguilar et al. (2018a), which favors the interpretation of the observed CR break around $300$ GV as a break in diffusion rather than CR injection:

[TABLE]

where $\delta$ and $\delta_{2}$ are, respectively, the indices below and above the break at $R_{1}$, $D_{0}$ is the overall normalization, and $\beta=v/c$ the velocity of the CRs. The velocity of convective winds, $\bm{V}(\bm{x})$, is assumed to be constant and orthogonal to the Galactic plane, $\bm{V}(\bm{x})={\rm sign}(z)\,v_{0,{\rm c}}\,{\bm{e}}_{z}$. Diffusive reacceleration is parametrized by the velocity $v_{\rm A}$ of Alfvèn magnetic waves Dogiel et al. (1990); Seo and Ptuskin (1994):

[TABLE]

This formula exhibits an ambiguity regarding the choice of $\delta$ when, as in our case, spatial diffusion has a break. We choose to use a single $\delta$ at all rigidities, specifically the $\delta$ below the break $R_{1}$ introduced in Eq. (5). Finally, CR propagation is affected by energy losses, continuous, adiabatic and catastrophic, for which we use the default Galprop implementation. Explicitly, in Eq. (1) continuous energy losses are described by the term $\mathrm{d}p/\mathrm{d}t$, and catastrophic losses through fragmentation or decay are described by the respective lifetimes $\tau_{f}$ and $\tau_{r}$. Catastrophic energy losses provide a source term for secondary CRs (see Sec. IV for details). Consequently, the propagation equations for the different CR species (primary and secondaries) constitute a coupled system of differential equations.

Beside secondary production by CR proton and helium interactions with the ISM, annihilation of DM in the Galaxy may lead to an additional source of antiprotons from the fragmentation and decay of the annihilation products Bergstrom et al. (1999); Donato et al. (2004); Bringmann and Salati (2007); Donato et al. (2009); Fornengo et al. (2014); Hooper et al. (2015); Pettorino et al. (2014); Boudaud et al. (2015); Cembranos et al. (2015); Cirelli et al. (2014); Bringmann et al. (2014); Giesen et al. (2015); Evoli et al. (2015, 2012); Reinert and Winkler (2018). The corresponding source term reads

[TABLE]

where $m_{\mathrm{DM}}$ and $\rho(\bm{x})$ are the DM mass and DM energy-density profile, respectively. The sum runs over all DM annihilation channels $f$. The corresponding velocity averaged cross section and antiproton energy spectrum per annihilation are denoted by $\left\langle\sigma v\right\rangle_{f}$ and $\mathrm{d}N^{f}_{\bar{p}}/\mathrm{d}T$, respectively, where $T$ is the kinetic energy. Note that the factor $1/2$ corresponds to Majorana fermion DM. The energy spectra per annihilation at production, ${\mathrm{d}N^{f}_{\bar{p}}}/{\mathrm{d}T}$, depend on the DM mass, the kinematics of the annihilation process and details of the fragmentation and decay of the annihilation products. In this article, we consider the annihilation into bottom quarks, ${\rm DM\;DM}\to b\bar{b}$, for illustration. Note that the antiproton spectra of other hadronic channels exhibit similar shapes. Therefore, we expect this choice to be suitable to analyze a potential DM signal even though the true composition of annihilation channels may be very different. The corresponding best-fit DM mass may, however, differ for different annihilation channels, as discussed in Cuoco et al. (2017b). We employ the spectra presented in Cirelli et al. (2011).

The DM density profile is still subject to sizable uncertainties. However, CR probes of DM do not exhibit a strong sensitivity on the chosen profile Cuoco et al. (2018). Within this work we use the Navarro-Frenk-White density profile Navarro et al. (1996), with a characteristic halo radius as $r_{\mathrm{h}}=20\,$kpc and normalization so as to obtain a local DM density $\rho_{\odot}=0.43\,$GeV/cm3 Salucci et al. (2010) at the solar position $r_{\odot}=8\,$kpc.

We assume that CRs are in a steady state, which is a good approximation for nuclei. The geometry of the Galaxy is approximated by a cylindrical box with an extension of $r=20$ kpc and $z=\pm z_{\text{h}}$. We solve the diffusion equation numerically with the Galprop code333http://galprop.stanford.edu/ Strong et al. (2000); Strong (2015). Galprop solves the equation on a grid in the kinetic energy per nucleon and in the two spatial dimensions $r$ and $z$. The radial and $z$ grid steps are chosen as $\Delta r=1\,$kpc and $\Delta z=0.1\,$kpc, respectively. The grid in kinetic energy per nucleon is logarithmic between $1$ and $10^{7}\,$MeV with a step factor of $1.4$. We use versions Galprop 56.0.2870 and Galtoollibs 855444https://galprop.stanford.edu/download.php as the basis and implement several custom modifications. Most importantly, first, we include the antiproton production cross sections from di Mauro et al. di Mauro et al. (2014) and Winkler Winkler (2017) (with updated parameters from Korsmeier et al. Korsmeier et al. (2018b)). We allow for either using default cross sections from tables or the full parametrization of the Lorentz invariant cross sections. In the latter case the Lorentz transformation and angular integration is performed on the fly. Secondly, we implemented the possibility of using a smoothly broken power law in the injection spectrum.

At low energies, CRs are deflected and decelerated by solar winds and the solar magnetic field. The strength varies in a 22-year cycle and is, therefore, commonly referred to as solar modulation. The effect starts to be significant below few tens of GV and increases towards low energies. Solar modulation can be described by a propagation equation similar to the one of interstellar CR propagation but adjusted to the situation in the heliosphere. It can be solved numerically and compared to data Kappl (2016); Maccione (2013); Vittino et al. (2018); Boschini et al. (2018). Progress in the understanding of CR propagation in the heliosphere has been made recently Tomassetti (2017); Cholis et al. (2016) especially thanks to (i) the data of the Voyager I probe which has left the heliosphere a few years ago and thus determines the CR fluxes before they are influenced by the solar effects, and (ii) the data of PAMELA, which provided time dependent CR fluxes. Further progress is expected from the analysis of the time-dependent CR fluxes released recently by AMS-02 Aguilar et al. (2018b, c). Nonetheless, a detailed understanding is still missing at the moment. For this reason we resort to the commonly used force-field approximation Fisk (1976), where the CR flux near Earth is calculated as

[TABLE]

Here, $e$ is the elementary charge, $Z_{i}$ and $m_{i}$ are the charge number and mass of the species $i$, respectively, and $\varphi_{{\rm SM},i}$ is the corresponding solar modulation potential. The variables $E_{\text{LIS},}$ and $E_{\oplus,i}$ describe the total CR energy before and after solar modulation, respectively. To minimize the deviations of solar modulation from the force-field approximation we limit the analysis to rigidities larger than 5 GV. In this respect various comments are in order:

• (i)

The recent publication by AMS-02 on time-dependent proton and helium fluxes Aguilar et al. (2018b) indicates that the proton to helium ratio is in general constant, except for data taken after May 2015 at rigidities below 3 GV. We will use proton and helium data before May 2015, and we can thus use the same modulation for the two species.

• (ii)

It is well established that solar modulation is charge-sign dependent. The AMS-02 measurements of the time-dependent proton and helium fluxes Aguilar et al. (2018b) as well as the electron and positron fluxes Aguilar et al. (2018c) provide a further confirmation of this effect. To take this effect into account we will use a different solar modulation potential for antiprotons as compared to $p$ and He, see Sec. III.

• (iii)

The AMS-02 time-dependent data show explicitly a deviation from the force-field approximation. For example, the proton spectra from different periods exhibit crossings,555For instance, the proton fluxes from November 2013 (Bartels rotation 2460) and March 2015 (Bartels rotation 2476) cross at $\sim 4$ GV. an effect that is not possible in the force-field approximation. This in particular is seen at low energies below $\sim 5$ GV and during the maximum of solar activity.

We conclude that below 5 GV the force-field approximation starts to lose its reliability. Hence we discard data below 5 GV by default in our analysis. Note also that the data we use do not include the solar maximum. Nonetheless, we also performed test fits with different lower cuts on the rigidity, still within the force-field approximation. From these fits we will further show that a cut of 5 GV is conservative from the point of view of DM searches.

III Default setup

Here we provide a short summary of the fit setup which, in general, is rather similar to the one used in Korsmeier and Cuoco (2016); Cuoco et al. (2017a). However, there are some important differences which we point out below. We use AMS-02 proton and helium fluxes Aguilar et al. (2015a, b), which both span the data period from May 2011 to November 2013, and the AMS-02 antiproton to proton ratio Aguilar et al. (2016), taken during the period May 2011 to May 2015. Furthermore, we use proton and helium data from Voyager Stone et al. (2013) and, in some fits, complement with data from CREAM Yoon et al. (2011).

In general, the likelihood for the CR fit is given by the product of the likelihoods of all experiments and CR species:

[TABLE]

Here $\phi^{(e)}_{{s},i}$ is the flux measured by the experiment $e$ for the CR species $s$ at the rigidity $R_{i}$ and $\phi^{(\text{m})}_{s,e}$ is the corresponding Galprop model. The covariance matrix $\mathcal{V}^{\left(e,s\right)}$ describes the uncertainty of the flux measurement. In the default setup we assume uncorrelated uncertainties, $\mathcal{V}_{ij}^{\left(e,s\right)}=\delta_{ij}\left[\sigma^{(e)}_{{s},i}\right]^{2}$. We suppress the explicit dependence of $\phi^{\text{(m)}}$ on all the fit parameters.

Cosmic-ray propagation is described by a total of 15 (or 17, when including DM) parameters. They are partly described in Sec. II, but, for convenience, we list them again below. Six parameters are used to describe the injection spectrum of protons and helium, i.e., the slopes below and above the rigidity break, $\gamma_{1,p}$,$\gamma_{1}$, $\gamma_{2,p}$, $\gamma_{2}$, the common rigidity break $R_{0}$ and a common smoothing parameter $s$. Five more parameters describe propagation, i.e., the normalization $D_{0}$ and slope $\delta$ of the diffusion coefficient, the velocity of Alfvèn magnetic waves, $v_{A}$, the convection velocity, $v_{0c}$, and the Galaxy’s half-height, $z_{\mathrm{h}}$. In the default setup we limit the fit range of AMS-02 data from 5 to 300 GV. A more detailed discussion and justification of these numbers is given further below. We notice, however, that in this way we avoid fitting the two parameters describing diffusion above 300 GV, which are instead fixed to $R_{1}=300$ GV and $\delta_{2}=\delta-0.12$. Two further parameters, $m_{\rm DM}$ and $\langle\sigma v\rangle$, are used to parametrize the considered DM model when antiprotons from DM annihilation are included in the fit. These 11 (13) parameters are scanned using MultiNest Feroz et al. (2009). For the MultiNest setup we use 500 live points, an enlargement factor efr=0.7, and a stopping criterion of tol=0.1. The final efficiency of the scan is typically around 7% and the number of likelihood evaluation around 200 000.

The remaining four parameters are treated in a simplified way. Two are the normalization of the proton and helium fluxes, $A_{p}$ and $A_{\mathrm{He}}$, respectively, and the other two are the solar modulation potential, $\varphi_{\text{SM,AMS-02},p,{\rm He}}$, referring to $p$ and He as well as $\varphi_{\text{SM,AMS-02},{\bar{p}}}$ referring to ${\bar{p}}$. We do not apply priors on $\varphi_{\text{SM,AMS-02},p,{\rm He}}$, while for $\varphi_{\text{SM,AMS-02},{\bar{p}}}$ we apply a very weak Gaussian prior, i.e., we add to the main likelihood the term $-2\,\log({{\cal L}_{{\rm SM}}})=(\varphi_{\text{SM,AMS-02},p,{\rm He}}-\varphi_{\text{SM,AMS-02},{\bar{p}}})^{2}/\sigma_{\varphi}^{2}$ where $\sigma_{\varphi}=100$ MV. We profile over these four parameters on-the-fly at each MultiNest likelihood evaluation following Rolke et al. (2005). More precisely, for each evaluation in the fit within the 11- (13-) dimensional parameter space the likelihood is maximized over the four remaining parameters. This maximization is performed with Minuit James and Roos (1975). To interpret the scan result we use a frequentist framework, and we build one- and two-dimensional profile likelihoods in the different parameters, from which we derive contours which are shown in various figures in the following.

There is a subtlety in the treatment of solar modulation. As mentioned above, the AMS-02 ${\bar{p}}/p$ ratio data – and hence the proton data used in this ratio – are taken from a different period than the AMS-02 $p$ and He data. Thus the two $p$ datasets should be modulated by a different solar modulation potential to be self-consistent. An obvious improvement would be to use $p$ and ${\bar{p}}/p$ datasets from the same time period, which are, however, not available in the AMS-02 publication Aguilar et al. (2016). Alternatively, the ${\bar{p}}$ absolute data (which are available in Aguilar et al. (2016)) could be fitted instead of the ${\bar{p}}/p$ ratio. This procedure is, however, suboptimal, because in the ratio some systematic uncertainty cancels out and thus the ${\bar{p}}$ data have a larger relative error than the ${\bar{p}}/p$ ratio data. Moreover, the use of the ${\bar{p}}/p$ ratio data considerably simplifies the fit since the ratio is considerably less sensitive to the injection parameters than the absolute ${\bar{p}}$ flux. To resolve this issue we proceed as follows. We derive an “effective” $p$ spectrum from ${\bar{p}}/p$ ratio and ${\bar{p}}$ flux from Aguilar et al. (2016) taking ${\bar{p}}/({\bar{p}}/p)$. We then divide the published $p$ spectrum Aguilar et al. (2015a) by this effective $p$ spectrum. As expected, this ratio is consistent with 1 above $\sim 40$ GV and slowly decreases at low rigidities due to solar modulation, up to a maximum deviation of 5% at 1 GV. As a function of rigidity this ratio is then approximated by a smooth log-parabola, which in turn is used to multiply the published $p$ and He data Aguilar et al. (2015a, b) to create an effective $p$ and He data set that corresponds to the same period as the ${\bar{p}}/p$ ratio. During our study the AMS-02 publication Aguilar et al. (2018b) became available providing time-dependent $p$ and He fluxes. This provides us with an alternative possibility to derive the $p$ and He effective fluxes corresponding to the period of the ${\bar{p}}/p$ ratio, namely, averaging the given monthly fluxes over that period. This second method is expected to be more robust. Nonetheless, we found that the effective $p$ and He fluxes build from the two methods perfectly agree, except below 3 GV where we find a difference of the order of 2%, which is comparable with the error bars. Using these effective $p$ and He fluxes we can self-consistently use the modulation potential $\varphi_{\text{SM,AMS-02},p,{\rm He}}$ for the absolute fluxes of $p$ and He and for the $p$ in the ${\bar{p}}/p$ ratio, and the modulation potential $\varphi_{\text{SM,AMS-02},{\bar{p}}}$ for the ${\bar{p}}$ in the ${\bar{p}}/p$ ratio. In the future, to avoid these issues, we recommend experimental collaborations to periodically release global datasets including the measurements of the fluxes of all species available and all referring to the same time period.

As a default production cross section for the secondary antiprotons we use the model from Winkler et al. Kappl and Winkler (2014); Winkler (2017) with updated parameters from Korsmeier et al. Korsmeier et al. (2018b) (referred to as param. MW in the following). The results of an alternative cross-section parametrization by di Mauro et al. di Mauro et al. (2014) (param. MD) and the effect of uncertainties of the cross-section parameters are discussed in detail in Sec. IV.

The results of the fit with this default setup, with and without DM, are shown in Figs. 1–3 and listed in the first two columns of Tab. 3. Figure 1 contains the comparison of the best-fit proton and helium spectra to data from AMS-02 and Voyager. The respective plots from the fit including DM look very similar. The residuals in the lower panels show a perfect agreement of the AMS-02 data with the Galprop model. However, they also already hint at a problem: The fluctuation of data points around the best fit is much smaller that the dominating systematic uncertainty. If this uncertainty is taken to be uncorrelated the fit results in a $\chi^{2}$ per degree of freedom (dof) much smaller than 1. We elaborate on possible correlation scenarios in more detail in Sec. V.

Figure 2 shows the best fit of the antiproton-to-proton ratio. Considering the fit regime from 5 to 300 GV there is a clearly visible improvement in the residuals, if DM is included in the fit. The significance in terms of $\chi^{2}$ difference between the fits excluding and including DM is $\Delta\chi^{2}\simeq 12.7$ which formally corresponds to $3.1\sigma$. For the considered $b\bar{b}$ channel an annihilation cross section around $10^{-26}\,\text{cm}^{3}/s$ and a mass around 75 GeV provides the best fit. As already found in Cuoco et al. (2017a) we point out that the DM signature constitutes a spectral shape which is very different from the astrophysical secondary or tertiary components.

Finally, we summarize the best-fit CR and DM parameters in the triangle plot displayed in Fig. 3. Note that, compared to previous results in Cuoco et al. (2017a) and due to the change of the standard cross section, the value of $\delta$ has increased to 0.42 and 0.38 in the case without and with DM, respectively. These values are in much better agreement with expectations from B/C data which points to 0.4–0.45 Genolini et al. (2015); Kappl et al. (2015). The additional, embedded plots display the profiles for solar modulation. There is a small but nonsignificant preference for a slightly smaller solar modulation potential for antiprotons. From a theoretical point of view it is not exactly clear whether we expect a larger or smaller solar modulation potential for protons or antiprotons. The behavior is expected to depend on the polarity of the solar magnetic field. Since there is a change of polarity in 2013, which is in the middle of the period of the AMS-02 measurement of the antiproton-to-proton ratio, the exact situation is unclear. Future more detailed analyses of the recent monthly data by AMS-02 will provide a better understanding of this issue.

In the following we will discuss step by step the difference in our default setup with respect to the previous analysis performed in Cuoco et al. (2017a). In particular there are three main points: We

1. (i)

removed data from CREAM, as it provides a source of tension with high-rigidity data from AMS-02 that is likely to stem from a systematic uncertainty in the normalization which were not properly addressed, 2. (ii)

changed the ${\bar{p}}$ production cross section from Tan&Ng Tan and Ng (1983) to the more recent param. MW, 3. (iii)

introduced a separate solar modulation potential for $p$ and He, with respect to ${\bar{p}}$ as described above.

Besides the change of the best-fit propagation parameters (in particular $\delta$) mentioned above, also the significance of the potential DM signal changes. In Cuoco et al. (2017a) $\Delta\chi^{2}$ was found to be at the level of 25, corresponding to a significance of $\sim 4.5\sigma$. We performed several additional fits to trace the origin of the difference with respect to the current analysis. The $\Delta\chi^{2}$ changes

1. (i)

from 25 to 21, when removing CREAM data and fitting data only below 300 GV, 2. (ii)

from 21 to 11, when changing the ${\bar{p}}$ cross section to param. MW, 3. (iii)

and finally from 11 to 12.7, when introducing separate modulation potentials.

This investigation shows that the cross-section parametrization has a potentially large impact on the DM significance. We will investigate the robustness of the signal against cross-section uncertainties in the next section.

In order to examine the compatibility of our new setup with CREAM data we perform one further fit. We include data from CREAM and allow for four additional fit parameters. First, we include the break position of the diffusion coefficient, $R_{1}$, and slope above this break, $\delta_{2}$, which are sampled by MultiNest. Secondly, we introduce two normalization parameters of the CREAM proton and helium data which we leave free in order to take into account a possible systematic uncertainty (especially in the energy scale, which is degenerate with an uncertainty in the normalization for a pure power law the spectrum, as it is the case for the CREAM data in the limited energy range $\sim 1$ TeV–10 TeV). These nuisance parameters relieve the above-mentioned tension present when the data are used with the nominal normalization. Going from the default fit to this extended setup with a total of 19 (21) parameters does not introduce significant complications since the two extra MultiNest parameters $R_{1}$ and $\delta_{2}$ are mainly determined by the data above 300 GV and are thus largely uncorrelated with the rest of the parameter space, while the two normalization CREAM parameters are profiled over similarly as described above. The result of the fit is included in Tab. 3 as the third and fourth columns. It justifies our choice $R_{1}=300$ GV and $\delta_{2}=\delta-0.12$ made above. The DM significance is similar to the default setup and only increases slightly to $\Delta\chi^{2}=15.1$, which corresponds to $3.5\sigma$.

Finally, as announced in Sec. II we briefly discuss the effect of taking into account data below 5 GV while still using the force-field approach to describe solar modulation. To this end we performed three additional fits following the above default setting (17 parameters) but with a rigidity cut at 3, 2 and 1 GV. The corresponding values for the $\chi^{2}$ are shown in Tab. 1. We observe a significant gradual worsening of the overall fit quality with decreasing rigidity, in particular between 2 and 1 GV. This applies for both fits with and without DM. It indicates that the model becomes less and less able to explain the data at lower energies, as expected from various effects like deviations from the force-field approximation, or, possibly, deviations from the simple scenario of convection-reacceleration. The significance of a DM signal is maximal for a energy cut of 3 GV with a $\Delta\chi^{2}$ of around 19 and then decreases to $\Delta\chi^{2}\simeq 11$ for a cut at 2 GV and $\Delta\chi^{2}\simeq 0$ for a cut at 1 GV, where, however, the overall fit quality is significantly worse as mentioned above. The cut at 5 GV, hence, does not maximize the significance of a possible DM signal. Note that the values of the other CR parameters do not change significantly with the different rigidity cuts, while the errors have a slight improvement. Thus, the estimation of the CR parameters seems robust with respect to the variation of the lower rigidity cut.

IV Antiproton cross sections

IV.1 Introduction

As we have seen in the previous section, the ${\bar{p}}$ production cross section has an important impact on the fit and on the significance of the DM signal. In this section we thus go more into depth in the investigation of this issue. The CR antiprotons in our Galaxy are dominantly produced by the interaction of CR protons and helium with the ISM in the Galactic disk. The source term for the CR projectile nucleus $i$ and the ISM nuclei component $j$ is given by

[TABLE]

Here $\phi_{i}$ is the CR flux, $n_{\mathrm{ISM},j}$ is the density of the ISM, and $\mathrm{d}\sigma_{i,j}/\mathrm{d}T_{\bar{p}}$ is the energy-differential cross section for antiproton production. The parameters $T_{i}$ and $T_{\bar{p}}$ denote the kinetic energy of the CR projectile and antiproton, respectively. High-energy experiments measure the fully differential cross sections which are usually stated in the Lorentz-invariant form, $E_{\bar{p}}\,\mathrm{d}\sigma/\mathrm{d}p_{\bar{p}}^{3}$. There are two different strategies to extract the energy-differential cross section from the experimental data. On the one hand, Monte Carlo generators are tuned to the data and afterwards used to extract the required cross section Kachelriess et al. (2015); Lin et al. (2017). However, at the moment these approaches lack consistency with data at either low or high energies, depending on the specific generator Donato et al. (2017). On the other hand, an analytic parametrization of the Lorentz-invariant cross section is fitted to the experimental data. In this case, the energy-differential cross section is obtained by, first, performing a Lorentz transformation to the frame where the ISM component is at rest and, secondly, an angular integration di Mauro et al. (2014); Kappl and Winkler (2014); Winkler (2017); Korsmeier et al. (2018b). This approach works reasonably well throughout the whole energy range of AMS-02, namely, from a rigidity of 1 to 400 GV. Therefore, we rely on the parametrization approach in the following. More details are given in Donato et al. (2017); Korsmeier et al. (2018b). We exploit the two parametrizations, param. MW Winkler (2017) (used in the default setup above) and param. MD di Mauro et al. (2014), for which we use the parameters updated to the most recent data from NA61 and LHCb as presented in Korsmeier et al. (2018b). The uncertainty on the antiproton source term, solely due to cross sections, is at about 5% above $T_{\bar{p}}=5$ GeV and increases to 10% below. To take this uncertainty properly into account we apply two different methods.

IV.2 Covariance matrix method

In the first method, we propagate the error of the cross-section parametrization to the flux of the CR antiprotons. This method was already suggested and applied in Ref. Reinert and Winkler (2018). The procedure works as follows. We use the covariance matrix of the cross-section fits from Ref. Korsmeier et al. (2018b) and sample $N=1000$ random parameter combinations $k$ and the corresponding antiproton source terms $q^{(k)}_{\bar{p}}$. From these we determine the covariance matrix of the relative source term $q^{(k)}_{\bar{p}}(R_{i})/q^{\text{(best\,fit)}}_{\bar{p}}(R_{i})$ at rigidities $R_{i}$ of AMS-02 data points $i$. It is given by:

[TABLE]

(In formulas, tables and figures we abbreviate cross section with XS.) We assume that the covariance matrices of the relative source term and of the relative flux are identical.666This is a good approximation since the relative uncertainty is invariant under propagation which is described by a linear differential equation. In other words, the covariance matrix of the antiproton flux is given by

[TABLE]

Here $\phi^{(\text{AMS-02})}_{{\bar{p}},i}$ is the antiproton flux measured by AMS-02 at the rigidity $R_{i}$. Accordingly, for the log-likelihood of the antiproton data in our fit the covariance matrix is replaced by:

[TABLE]

Obviously, this method has one weakness: Error propagation in terms of a covariance matrix assumes that the likelihood of the original cross-section parameters and the corresponding likelihood in the space of AMS-02 antiproton flux data points are well approximated by a multivariate Gaussian distribution. However, the true likelihood might be more complicated. Furthermore, the method assumes that there is no correlation between the cross-section uncertainties and CR propagation uncertainties. To take into account these shortcomings we thus consider a second method below.

IV.3 Joint fit method

In the second method we perform a joint fit of CRs and the antiproton production cross section. By simultaneously fitting the CR and cross-section parameters to CR fluxes and experimental cross-section data, the full likelihood is correctly taken into account. The price that we have to pay is an increase in the number of parameters for an already high-dimensional and hence time-consuming fit. In our default setup the CR fit contains 11–13 MultiNest parameters; the fit of cross section uses 7–10 free parameters. Thus, naively merging the fit would lead to ${\cal O}(20)$ free MultiNest parameters, which is extremely challenging. Therefore, we reduce the number of free parameters in the cross-section parametrization to those which affect the shape of the CR antiproton source term the most (as discussed further below) and fix the remaining parameters. We focus on param. MW of the Lorentz-invariant cross section since the meaning of the single parameters is more obvious compared to param. MD. The parametrization MW for prompt antiproton production in proton-proton collisions depends on the center-of-mass energy $\sqrt{s}$ of the collision, the energy of the antiproton divided by the maximal antiproton energy $x_{\mathrm{R}}$, the transverse momentum of the antiproton $p_{\mathrm{T}}$, and the six parameters $\mathcal{C}=\{C_{1}...C_{6}\}$:

[TABLE]

where $m_{\text{T}}=\sqrt{p_{\text{T}}^{2}+m_{p}^{2}}$. The inelastic cross section $\sigma_{\mathrm{in}}$ of $pp$ scattering is defined in Winkler (2017). The factor

[TABLE]

describes the scaling violation of the cross section at low $\sqrt{s}$, and $X$ is defined by

[TABLE]

For nonproton nuclei in the projectile CR or target ISM state we rescale the $pp$ cross section as described in Korsmeier et al. (2018b). Furthermore, the total antiproton source term includes antiprotons produced by the decay of intermediate antineutrons or antihyperons. We apply the scalings from Winkler (2017). The total likelihood for the joint fit is given by the product of the CR and cross-section likelihoods:

[TABLE]

The procedure to fit the cross-section data follows Korsmeier et al. (2018b). We fit to the same datasets (NA49 Anticic et al. (2010), NA61 Aduszkiewicz et al. (2017), Dekkers et al. Dekkers et al. (1965), NA49 ($p$C) Baatar et al. (2013), LHCb ($p$He) Graziani et al. (2017)777During our analysis LHCb published the final analysis Aaij et al. (2018) of the cross section. They differ from the preliminary results by a scale factor of about 10%. However, since we include a scale uncertainty of 10% in our analysis, we do not expect a significant effect on the results.) and use the same likelihood definition:

[TABLE]

As before, $e$ denotes experiments (this time of cross-section measurements) with data points $i$, while $m$ denotes the cross-section parametrization. The symbols $\sigma_{\text{inv}}$ and $\sigma_{\sigma_{\text{inv}}}$ represent the Lorentz-invariant cross section and its uncertainty, respectively. We account for a scale uncertainty $\omega_{e}$ of each cross-section measurement $e$. It is constrained by a Gaussian prior, namely, the second term in Eq. (19). During the fit, these parameters are treated in a simplified way, similar to the normalizations and solar modulation parameters (cf. Sec. III). For each MultiNest step the $\omega_{e}$s are profiled over by performing a Minuit fit.

We now discuss our choice of parameters considered in the fit. The uncertainty of the antiproton production cross section has different origins. At high energies, above $\sqrt{s}=10$ GeV, the shape of the cross-section data is constrained extremely well. The largest uncertainty is the normalization of the cross section. The origin of this uncertainty is the experimental difficulty to determine the luminosity better than a few percent. At lower energies, data is more scarce and less precise. Furthermore, the theoretically motivated and experimentally confirmed concept of scaling invariance of the cross section is broken. Therefore, extrapolations are less trustworthy. We thus identify the parameters $C_{1}$ and $C_{5}$, $C_{6}$ as the most relevant for our purpose. They determine the normalization of the whole cross-section parametrization and the shape at low energies, respectively. In the following joint fit, we vary those parameters only, while all the other parameters are chosen to be fixed to the values from Korsmeier et al. (2018b). Note that adding these three parameters and the cross-section data to the global fit complicates the structure of the likelihood. Using the same configuration of MultiNest the total number of likelihood evaluations increases up to 500 000. This is a clear disadvantage with respect to the approach with a covariance matrix which, instead, does not significantly complicate the fit.

IV.4 Results

In total, we perform eight fits to study the effect of cross sections and their uncertainty on the results. These corresponds to four different setups with respect to the cross-section parametrization, while for each setup we perform one fit with and one without DM. The results of all fits are summarized in Tab. 4. The residuals of the antiproton spectra are shown in Fig. 4. We start by discussing the effect on the fits without DM. In the first setup, we change to cross-section param. MD. As a result the fit quality improves marginally from $\chi^{2}_{\rm min}=35.6$ to 34.2. The best-fit parameters are compatible and their uncertainty is very similar. In Fig. 5 we show the comparison of the fit contours with our default setup. The contours of this setup are slightly increased, in particular towards a larger $v_{0,{\rm c}}$, but they are not systematically shifted. Then, we apply the covariance matrix method to both cross-section parametrizations. We find that the effect on the best-fit parameters and their uncertainties is negligible. It is, however, interesting to observe that in the residual plot of this fit there is a systematic shift of all points towards larger values. This is because in the energy range above 5 GV the cross-section covariance matrix $\mathcal{V}_{\text{XS}}^{\left(\phi^{\text{AMS-02}}_{\bar{p}}\right)}$ mostly encodes the cross-section normalization uncertainty.

The results of the joint fit method are compatible with both the previous two fits, the one without cross-section uncertainties and the covariance matrix method. The residuals are mostly unchanged, and the parameter estimation and the corresponding uncertainty are similar. Nevertheless, the joint fit provides very interesting insights allowing for further cross-checks. In particular, it allows us to investigate whether the CR parameters and the cross-section parametrization are correlated. In the extreme case, one might even imagine that the very precise CR data could constrain the cross-section parametrization. Figure 6 displays the part of the parameter triangle which shows the correlation of the cross-section parameters with all CR parameters. We conclude that there is no significant correlation between the CR and cross-section fit parameter. Consequently, we expect that the cross-section parametrization is not affected by the CR data. This is confirmed in Fig. 7, which shows the $\Delta\chi^{2}$-profiles of the three cross-section parameters varied in the fit. We compare the profile of the total likelihood with the profile of a fit to only the cross-section data: Both profiles agree well within their respective uncertainties. In fact, only $C_{5}$ is shifted to slightly lower values in the joint fit. This absence of correlation likely explains also why the covariance matrix method performs reasonably well and gives similar results.

The impact of the cross-section uncertainty on the possible DM hint in the antiproton spectrum can be understood looking at the residuals in Fig. 4 for the case of the fits with DM. In all scenarios the flattening of the residuals is similar. In terms of $\chi^{2}$s the improvement of the fit with DM compared to the fit without DM was $\Delta\chi^{2}=12.7$ for the case of a fixed MW parametrization (default setup). The covariance matrix and joint fit methods decrease the $\Delta\chi^{2}$ to 10.9 and 10.7, respectively, indicating that the evidence for DM is not strongly affected by the cross-section uncertainties. Furthermore, the result of the best-fit DM mass and velocity-averaged annihilation cross section is not strongly affected by the uncertainties. We show the comparison of the best-fit contours in Fig. 8. For comparison we also show the cross-section limit derived from gamma-ray observations of dwarf spheroidal galaxies Albert et al. (2017) and the best-fit region of the GCE Calore et al. (2015) for the considered $b\bar{b}$ channel. All observations are compatible, in particular, since they are affected by astrophysical uncertainties in different ways providing additional freedom to alleviate a certain tension among them, see Cuoco et al. (2017b) for a detailed analysis of the subject.

Above, we have focused on fits and results where we exclude data below 5 GeV, since as argued in the introduction, the results using data down to 1 GV are more prone to further systematic uncertainties, especially solar modulation. Nonetheless, it is interesting to have a look at the fit results including the low-energy data from a methodological point of view. The cross-section (shape) uncertainties are most severe at low energies, while at higher energies only the normalization is uncertain. Therefore, it is not very surprising that the results of both methods are very similar. If, however, we include data at low energies, the picture changes. We investigated how the best-fit parameters are affected by the two methods and find that both methods have a still small, but similar effect on the parameter space. Furthermore, we observe that the error contours of the covariance matrix method are a bit larger compared to the joint fit method, in other words, the former is more conservative. We regard this as proof of concept: The covariance matrix method, which is easier to implement and less time consuming in the fit, is a reasonable approximation to the more complete joint fit method.

The above results are somehow at odds with the results of Reinert and Winkler (2018), where flat ${\bar{p}}$ residuals are achieved down to 1 GV and no significant preference for a DM signal (a global significance of $1.1\sigma$) was found. The authors of this study use a covariance matrix method to account for the cross-section uncertainties. They conclude that the inclusion of these uncertainties is the main reason why their analysis does not provide a hint for DM. Nonetheless, the results shown above indicate that the cross-section uncertainties do not have such a strong impact. An important difference is that in Reinert and Winkler (2018) only the ${\bar{p}}$ spectrum is fitted, with the source terms for ${\bar{p}}$ being fixed using the observed $p$ and He spectra corrected for solar modulation. This has the advantage that the injection parameters do not need to be fitted, although it requires some assumption on how to extrapolate the observed local $p$ and He spectra to the ones for the whole Galaxy needed for the secondary source terms. Instead, in our approach ${\bar{p}}$, $p$ and He are fitted simultaneously and we include $p$ and He injection parameters in the fit. Fitting the $p$ and He spectra provides extra constraints on the propagation with respect to fitting ${\bar{p}}$ only. For example, it is well known, e.g., Moskalenko et al. (2002), that strong reacceleration produces a low-energy ($\lesssim 10$ GeV) bump in the $p$ spectrum, which is not observed. The $p$ spectrum, thus, provides strong constraints on the amount of reacceleration, although this is, in part, degenerate with the break in the injection Moskalenko et al. (2002). We thus suspect that in Reinert and Winkler (2018) it is possible to accommodate the secondary ${\bar{p}}$ spectrum, while this is not possible anymore when constraints from $p$ and He are included as it is the case in our analysis. Further differences concern a different treatment of reacceleration (which in Reinert and Winkler (2018) is confined to the Galactic disk only, while it is uniform over the whole diffusion region in our case), adiabatic energy losses from convection and a two-dimensional source term distribution used in our analysis. Therefore a direct comparison is not easily achievable and would require a substantial modification of our setup, which is left for future work.

V AMS-02 correlations

With the era of space-based CR detectors the statistics and quality of collected data have significantly increased. This also means that the relative weight of systematic uncertainties with respect to the statistical error has become more important. For example, the error budget of the measured proton and helium spectra is now completely dominated by systematics in most of the energy range. The question of how to assess and treat these uncertainties in a statistically correct way has thus become more pressing. The commonly used strategy is to add statistic and systematic uncertainties in quadrature and consider the new errors uncorrelated in energy. This, as we have seen in the previous sections, typically results in quite low $\chi^{2}$ values of the fit to the data. Table 3 shows that the typical values of the $\chi^{2}/\text{dof}$ are of the order 0.06–0.1 for proton and helium and 0.25–0.50 for the antiproton-over-proton ratio, whereas for a consistent treatment of all uncertainties one would expect a $\chi^{2}/\text{dof}$ close to 1. We conclude that either the systematic uncertainties of AMS-02 are overestimated or that there is a sizable correlation of the systematic uncertainty which is not correctly taken into account. The fact that our $\chi^{2}$ values are significantly smaller than expected has potentially problematic consequences: First, we cannot say anything about the goodness of fit and, secondly, the uncertainties deduced for our fit parameters could be affected. Thirdly, the significance of the potential DM signal may depend on the correlations in the systematic uncertainties. A proper assessment of this problem would require the experimental collaborations to provide the covariance matrix of the data points based on the knowledge of the detector. Since this is not available at the moment, our goal below is to use a simple approach and gain an approximate understanding of the effects of neglecting correlations.

V.1 Methodology

We follow two different strategies to answer the above questions. In the first one, we simply set the systematic uncertainty to 1% of the flux or flux ratio before adding it to the statistical uncertainty in squares.888A similar approach was taken in Cavasonza et al. (2017). The typical systematic uncertainty stated in the AMS-02 publications of proton, helium, and antiproton fluxes is on the order of a few percent. The remaining systematic uncertainty is thus assumed to be fully correlated among the rigidity bins, i.e. equivalent to an overall normalization uncertainty. Since normalizations are already profiled over in our fit setup, this kind of uncertainty is already taken into account and does not need to be included again as an uncertainty in the data. Thus, in practice, this approach implies a significant reduction of the uncertainty with a potentially equally significant effect.

The above approach, however, does not address the question about the presence of shorter-range correlations among the data points. The most complete approach to the problem would be to start from the knowledge of the detector and model the systematic uncertainties which contribute to the total error budget. In the AMS-02 publications the various contributing systematics are listed. They contain the acceptance uncertainty, trigger uncertainty, rigidity scale uncertainty, and uncertainty from energy unfolding. By studying the single systematic effects and modeling them, it would be possible to build the covariance matrix of the data. Unfortunately, again, this requires an inside knowledge of the detector which is not publicly available. We thus resort, in the following, to a simpler, data-driven approach which is expected to model the effect of correlated uncertainties reasonably well. With this approach we aim at determining the systematic uncertainty and a possible correlation between data points. The focus of this strategy, in particular, is on the study of a possible short-range (in rigidity) correlation component. This component is potentially more critical for our analysis since it can affect the significance of sharp features like the ones expected from DM. On the contrary, as argued above, long-range correlations are basically equivalent to a normalization uncertainty and have a small impact. As a first step, we thus split the covariance matrix of the CR datasets into a sum of three parts:

[TABLE]

Here, the first part ${\cal V}_{\text{stat}}$ contains all statistical, i.e. uncorrelated, uncertainties. In other words, the entries of the correlation matrix for the $i$-th and $j$-th data point are given by ${\cal V}_{\text{stat},ij}=\delta_{ij}(\sigma_{\text{stat},i})^{2}$. The statistical error alone is generally available in the corresponding publications. The third part ${\cal V}_{\text{long}}$ describes the long-range correlations. Typical examples which would fall into this kind of uncertainty are normalization or tilts of the whole data set. The second part is ${\cal V}_{\text{short}}$ which describes the correlation of up to a few neighboring points. Our ansatz is

[TABLE]

where the three parameters ${\ell_{\text{corr}}}$, $f$, and $\alpha$ describe the correlation length (in terms of the distance in rigidity-bins),999The ansatz of Eq. (21) assumes that the correlation length does not depend on the rigidity bin. We notice that the rigidity binning chosen by the AMS-02 experiment is inversely proportional to the energy resolution of the instrument. So, for example, an uncertainty in the energy unfolding is expected to be described by this kind of the covariance matrix. the fraction of the systematic uncertainty which is correlated, and the shape of the correlation matrix, respectively. Our goal is to determine the three parameters from the data themselves. In practice, we regard the data as a realization of the true covariance matrix, and we try to reconstruct them assuming the above parametrization, using standard statistical inference. To this end, besides the covariance matrix, we also need a model of the true energy spectrum, which we take as a smooth multiply broken power law with three breaks for the antiproton-to-proton ratio and the helium flux, and four breaks for proton flux. The parametrization is similar to Eq. (3) but extends to a higher number of breaks. Then, in principle, our inference should proceed with a fit of the covariance matrix parameters and spectrum parameters together to the observed data. Here, instead, we apply a simplified two-step approach. In the first step, we fit the smoothly broken power law to the data using as $\chi^{2}$ the full systematic uncertainties assumed to be uncorrelated. In practice, this step is equivalent to absorbing the long-range correlated uncertainties into the smooth spectrum, and it is insensitive to the exact errors used to define the $\chi^{2}$. In the second step we thus fix the smooth energy spectrum to the one derived in the above step, and we use the residual with respect to the data points, ${\bm{x}}$, to constrain the covariance matrix. The log-likelihood for the parameters ${\ell_{\text{corr}}}$, $f$, and $\alpha$ is given by

[TABLE]

where the first term comes from the normalization of the multivariate Gaussian (in the data) which we use as likelihood. We checked with a toy Monte Carlo that the method is self-consistent. To this end, in the Monte Carlo we choose benchmarks values of ${\ell_{\text{corr}}}$, $f$, and $\alpha$, i.e., a benchmark covariance matrix, and from that we draw random values of ${\bm{x}}$. We then verify that from the above likelihood we correctly reconstruct the input values of ${\ell_{\text{corr}}}$, $f$, and $\alpha$ within uncertainties.

Figure 9 shows the results of this procedure applied to proton, helium, and antiproton-over-proton data. The figure shows the likelihood profiles and the corresponding 1$\sigma$ to 3$\sigma$ contours in the frequentist interpretation. As can be seen, the resulting constraints on the three parameters are not very strong. In particular, the correlation length ${\ell_{\text{corr}}}$ is not constrained at the 2$\sigma$ level for all three datasets, whereas $\alpha$ is constrained to small values for {\ell_{\text{corr}}}\raisebox{-3.99994pt}{ \overset{>}{\sim} }5 and unconstrained for smaller correlation lengths. The only well-constrained quantity is the fraction of the systematic uncertainty with a small correlation length, $f$. The fits converge to $f\sim 0.07$ and 0.1 for proton and helium, respectively, and f\raisebox{-3.99994pt}{ \overset{<}{\sim} }0.3 for antiprotons . We note that, interestingly, the potential correlation of all three datasets can be described by the same parameters within uncertainties. Experimentally, however, it is unclear if equal systematics properties are expected for the three cases. Therefore, we will consider the potential correlation independently for each CR species in the following.

Given the large uncertainty in the determination of ${\ell_{\text{corr}}}$, and $\alpha$, we choose three benchmark cases compatible with the constraints to study the effects on the fit results in more detail. Specifically, the benchmark scenarios correspond to three different fixed values for the correlation length, ${\ell_{\text{corr}}}=0,\,5$ and 10. Let us first consider proton and helium. Figure 9 shows that the values for $f$ and $\alpha$ are constrained at every fixed value of the correlation length. Hence, we fix their values to the maximum of the likelihood profiles at each ${\ell_{\text{corr}}}$. The exact values are summarized in Tab. 2. For the antiproton-to-proton ratio the situation is different for two reasons: First, the number of data points is smaller compared to the proton and helium datasets. Secondly, and probably even more importantly, the size of the statistical uncertainty is significantly larger. While proton and helium are clearly dominated by only systematic uncertainties, the size of systematic and statistical uncertainty is comparable for the ${\bar{p}}/p$ ratio. As a result, the best-fit $f$ converges towards 0 and is only constrained from above, while $\alpha$ is completely unconstrained. Since our goal is to study the effect of possible correlations in the AMS-02 data, choosing $f=0$, i.e., the maximum of the likelihood is not a good choice. Instead, to be conservative, we choose the maximum $f$ value allowed by our constraints. In particular, we choose $f=0.3$, which is approximately the maximum allowed value at a confidence level of 90%. For $\alpha$ we choose a representative value of 1. Again, these values are summarized in Tab. 2.

V.2 Results

In total, we perform eight fits to investigate the effect of correlations in the AMS-02 data. They correspond to four different setups of the uncertainties where each setup is fitted once with and without DM. Summarizing, the four setups are the fixed 1% systematic uncertainty without any correlation of data points, and the three different benchmark scenarios from the data-driven covariance matrix with correlation lengths 0, 5, and 10. The results of the various fits are shown in Tab. Scrutinizing the evidence for dark matter in cosmic-ray antiprotons. As for the case in which we studied the cross sections, we show the residuals of the antiproton-to-proton ratio in Fig. 10. The shown uncertainties are the square roots of the diagonal elements of the covariance matrix and, therefore, some care is needed when drawing conclusions directly from the figures in the cases of ${\ell_{\text{corr}}}=5$ or 10. In general, however, we can clearly observe that the uncertainties are considerably smaller with respect to the default fit. From a methodological point of view the smaller uncertainties result in a more complicated fit. In order to converge, the new fits require up to 1.8 million MultiNest likelihood evaluations using the same configuration as in the default setup. The fit with fixed systematic uncertainty at 1% converges to a $\chi^{2}$/dof of $77.4/143\approx 0.5$. In more detail, the $\chi^{2}$ of the AMS-02 antiproton-to-proton ratio is 44 which now matches well with the 42 ${\bar{p}}/p$ data points included in the fit, and it is thus reasonable also in terms of goodness of fit. On the other hand, the $\chi^{2}$’s of proton and helium are 8 and 10, respectively, and are still much smaller than the number of data points, i.e., 50 each. This hints to the fact that the uncertainty of antiprotons is estimated reasonably well while the one for proton and helium is still overestimated. The effect on the CR parameter estimation is shown in Fig. 11 (upper left panel): Compared to our default fit, the size of the uncertainty contours slightly shrinks. The only further notable effect is that the Alfvèn velocity is shifted to slightly larger values, but by less than 1$\sigma$. Adding DM in this scenario improves the fit by a $\Delta\chi^{2}$ of 30 which would correspond to an evidence of 5.1$\sigma$.

The three benchmark scenarios from the data-driven approach converge to a total $\chi^{2}$s in the range 250–280, which is about a factor of 2 larger than the dof of 143. The $\chi^{2}$ for ${\bar{p}}/p$ is in the range 40–70, which, again, is reasonable in terms of goodness of fit, especially considering the DM fits, where the $\chi^{2}$ range shrinks to 40–50. Instead, the main contribution to the total $\chi^{2}$ comes from $p$ and He which are in the range 80–100. With an $f\sim 0.1$ for $p$ and He the total error is at the level of $\sim$ 0.2%, so these high values of $\chi^{2}$ are probably not unreasonable since it is unlikely that the accuracy of the model is at such a high level of precision. Finally, from the triangle plots in Fig. 11 we see that the size of the uncertainty contours shrinks significantly with respect to the default fit, an outcome which is expected given the smaller error bars.

As already observed in the case of the first method, we notice a shift of $v_{\mathrm{A}}$ to slightly larger values although compatible at 1–2$\sigma$ level. Furthermore, there is now a preference for a small halo height which converges towards 2 kpc which is the boundary of the fit range. This small value of $z_{\mathrm{h}}$ is disfavored by observations of electron and positron spectra Di Mauro et al. (2014), but still possible. The result also shows that with very small errors the fit becomes sensitive to $z_{\mathrm{h}}$ even without the use of cosmic clock data like 10Be/9Be. This further stresses the importance of studying in detail the properties of the data to understand if such small uncertainties can be used and exploited. All three benchmarks show evidence for a DM signal. The significance depends on the correlation length, namely, the improvement in the fit quality reaches $\Delta\chi^{2}=34$ (corresponding to $5.5\sigma$) in the case of ${\ell_{\text{corr}}}=0$, while it drops to around 18 ($3.8\sigma$) for ${\ell_{\text{corr}}}=5$ and 10.101010We checked that these significances do not strongly depend on the value for $z_{\mathrm{h}}$, e.g. at a fixed value of $z_{\mathrm{h}}=4$ kpc we get a $\Delta\chi^{2}$ of 36, 20, and 21 for the cases of ${\ell_{\text{corr}}}=0$, 5, and 10, respectively. Again, the results with the above different setups and benchmark scenarios show the importance of modeling correlations correctly. Depending on the “true” covariance, results may be driven in different directions, as exemplarily pointed out here for the case of the potential DM signal. In particular, typically, we see that the use of covariance improves the significances and reduces the errors on the estimated parameters, which indicates that the full potential of the data is not yet exploited.

The shrinking of the contours with respect to the case in which the correlation among the data-points is neglected is also observed for the DM mass and velocity-averaged annihilation cross section. Figure 12 shows the contours of the benchmark scenarios of 1% fixed systematic uncertainties and ${\ell_{\text{corr}}}=5$, respectively. In the latter case, the contour is additionally shifted to slightly smaller values of the DM mass and to larger $\langle\sigma v\rangle$. The contours for ${\ell_{\text{corr}}}=0$ and 10 are very similar to ${\ell_{\text{corr}}}=5$. The fact that the contours do not extent to lower values of $\langle\sigma v\rangle$ as for the default fit case is directly related to the preference of the fit for small $z_{\mathrm{h}}$ values, which give a lower DM signal. Lower values of $\langle\sigma v\rangle$ in the default fit are, correspondingly, associated to larger values of $z_{\mathrm{h}}$ which provide an increased normalization of the DM signal. We remark, nonetheless, that the DM contours of all our fits are compatible with each other. We also stress that the contours are derived for a fixed value of the local DM density of $\rho_{\odot}=0.43$ GeV/cm3, as mentioned in Sec. II, and that the additional uncertainty in $\rho_{\odot}$ (conservatively in the range 0.2–0.7 GeV/cm3 Salucci et al. (2010)) would contribute to extend to contours to both smaller and larger values of $\langle\sigma v\rangle$.

Furthermore, we note that a combination of cross-section uncertainties and correlations in the AMS-02 data does not significantly change the picture. We checked this by performing a fit taking into account the uncertainties of the cross-section parametrization (following the joint fit method as introduced in Sec. IV.3) as well as taking into account the 1% uncorrelated uncertainty from AMS-02. As expected, the resulting fit is very similar to the one fixing the cross-section parameters to their best-fit values.

Finally, we remark again that all results of this section have to be taken with caution. By guessing different possible scenarios of the AMS-02 uncertainties and its covariances in data, we demonstrate, merely as a proof of concept, that correlations in data can impose an important effect on the results, in particular, on the significance of the potential DM signal. Nonetheless, we do not deem our four benchmark scenarios as fully representative of the complete set of possible covariances in data or to contain the true covariance matrices. Our main aim is, instead, to show the importance of these effects in order to trigger further investigations on the subject, especially from the experimental collaborations, which, through the knowledge of the detector, have handles to study this issue more in detail. As we have shown, a better characterization of the covariance of the uncertainties can provide much more constraining power, and thus fully exploit the potential of the experimental data.

VI Conclusion

In this paper we investigated two important sources of uncertainties in the interpretation of global CR fits to AMS-02 data: the antiproton production cross section and the correlations in the experimental data. In particular, we investigated their impact on the potential dark matter signal found in previous studies Cuoco et al. (2017a); Cui et al. (2017); Cuoco et al. (2017b).

First, we reproduced the findings of Cuoco et al. (2017a) in an updated setup using the most recent cross-section parametrization from Korsmeier et al. (2018b) and an improved treatment of solar modulation. In addition, we studied a possible tension between the AMS-02 and CREAM data at large rigidities by either excluding CREAM data from the global fit or by taking into account a possible normalization uncertainty. We found that the potential DM signal persists, although with a slightly smaller significance of around $3\sigma$ (with respect to $4.5\sigma$ found in Cuoco et al. (2017a)). The reason for the reduction in significance is mainly due to the updated antiproton production cross-section parametrization.

We studied the uncertainties induced by the antiproton cross-section parametrizations using two approaches. On the one hand, we described the uncertainties by a covariance matrix assuming a multivariate Gaussian distribution. On the other hand, we performed a joint fit of CR and the most relevant cross-section parameters. We found that both descriptions lead to comparable results. As the latter approach reveals, this is due to the absence of strong correlations between the CR and cross-section parameters. We hence concluded that neither the results for the CR propagation parameters nor the significance of the DM signal are strongly affected by the antiproton cross-section uncertainties.

In all fits the $\chi^{2}/\text{dof}$ is significantly below one. We expect this to be a result of neglecting correlations in the AMS-02 data, which are not provided by the experimental collaboration. A long-distance correlation (in rigidity bins) amounts to an overall normalization or tilt of the spectrum, and it can be absorbed into the fit parameters. Using the systematic errors as given by AMS-02 and assuming these to be uncorrelated can, hence, greatly overestimate the error. We followed different approaches in order to illustrate the importance of the knowledge of these correlations. On the one hand, we simply reduced the uncorrelated error to 1% assuming the remaining error to be fully correlated (and absorbed in the overall normalization). As a result the overall fit quality worsens slightly, however, still providing $\chi^{2}/\text{dof}$ well below one. Interestingly, the preference for DM becomes more significant reaching $5\sigma$. On the other hand, we followed a data-driven approach, determining the maximal fraction of a correlated uncertainty for three choices of the correlation length (in terms of the distance between rigidity bins), $\ell_{\text{corr}}=0,\,5,\,10$. Our analysis shows that only a small fraction of the systematic uncertainty can be correlated this way, namely, about 10% for the $p$ and He data and no more than 30% (at 90% C.L.) for the ${\bar{p}}/p$ data. The remainder of the systematic uncertainty is compatible with just an overall shift of the global normalization. However, the identification of the exact shape of a covariance matrix from data and a typical correlation length is not unique. As an illustration and as a proof of concept, we choose three benchmark scenarios, although these are not expected to cover the whole ranges of possibilities. A more realistic characterization of the uncertainties can only be achieved in the future with further insight from the AMS-02 detector and analyses. Since the fraction of short-range correlated uncertainties is very small, for the moment treating AMS-02 systematic uncertainties as uncorrelated (i.e. statistical) is conservative. Nonetheless, we have shown that with a proper treatment of the uncertainties, the data can be much more constraining, both on the determination of propagation parameters and on the presence of a signal for DM. In particular, we find that the significance of a DM annihilation signals increases up to 5.5$\sigma$ depending on how we model the correlation. We thus stress that getting this kind of insight in the future would be extremely valuable and would allow one to better exploit the potential of the experimental data.

Finally, we note that two aspects remain to be studied in more detail in the future. First, solar modulation is not fully understood for rigidities below about $5$ GV and has to be investigated using the recently published monthly data of AMS-02. In addition, different diffusion models have to be studied to allow to include other secondary CR data such as 3He, Be, Li and B into our global analysis.

Acknowledgments

We thank Henning Gast, Björn Sarrazin and Stefan Schael for helpful discussions and comments. J.H. acknowledges support from the F.R.S.-FNRS, of which he is a postdoctoral researcher. The work of M. Korsmeier is supported by the “Departments of Excellence 2018-2022” grant awarded by the Italian Ministry of Education, University and Research (MIUR) (L. 232/2016). Simulations were performed with computing resources granted by RWTH Aachen University under projects rwth0085 and rwth0272.

Note added: Shortly after this analysis was made public, Refs. Cholis et al. (2019); Lin et al. (2019) appeared on the arXiv discussing similar subjects. The authors of Cholis et al. (2019) investigate the impact of cross-section uncertainties as well as time, charge and energy-dependent solar modulation effects on the significance of the DM signal, and arrive at conclusions which are similar to ours. In Lin et al. (2019) the authors exploit machine learning to marginalize more efficiently over cosmic-ray propagation uncertainties.