Fast and accurate AMS-02 antiproton likelihoods for global dark matter fits

The term $q_{i}(\bm{x},p)$ in eq. (2) incorporates the injection of CRs at each Galactic position $\bm{x}$, with a given energy dependence. Depending on the considered CR species $i$, we include primary and/or secondary contributions, which are briefly illustrated in the next subsections. We define as primary CR sources both the standard astrophysical sources such as supernova remnants (section 2.1.1), along with more exotic CR injections such as the production of (anti)particle CRs in the annihilation or decay of DM particles in the Galactic DM halo (section 2.1.2). Secondary CRs are instead produced by the interaction of the primary CRs with the gas in the interstellar medium (ISM) through fragmentation or decay processes (section 2.1.3).

As outlined above, the source term of CR nuclei significantly differs between primary and secondary CRs as well as between DM and astrophysical sources. In contrast, the mechanism and, therefore, the description of CR propagation are the same for all species. The dominant process for CR nuclei, especially at high energies, is the scattering on the turbulent magnetic fields in our Galaxy. Effectively, this leads to diffusion of CRs in a halo extending a few kpc above and below the Galactic plane. From linear perturbation theory, the diffusion coefficient is expected to be antiproportional to the spectrum of magnetic wave turbulence [75]. The spectrum of magnetic turbulence depends on the exact turbulence model; the typical benchmarks are Kolmogorov or Kraichnan, both leading to power-law dependence of the diffusion coefficient as a function of rigidity, but with different spectral indices of $\delta=0.33$ and 0.5 [76], respectively. We use a phenomenological approach and model the diffusion coefficient as a double-broken power law in rigidity

$\displaystyle D_{xx}\sim\beta R^{\delta_{l}}\cdot\left(1+\left(\frac{R}{R_{D,0}}\right)^{1/s_{D,0}}\right)^{s_{D,0}\,(\delta-\delta_{l})}\cdot\left(1+\left(\frac{R}{R_{D,1}}\right)^{1/s_{D,1}}\right)^{s_{D,1}\,(\delta_{h}-\delta)}\,,$

(2.5)

which is normalized to $D_{xx}(R=4\;\mathrm{GV})=D_{0}$. The spectral indices below, between, and above the two breaks at $R_{D,0}$ and $R_{D,1}$ are labeled $\delta_{l}$, $\delta$, and $\delta_{h}$. At the positions of the breaks, the power law is smoothed by the parameters $s_{D,0}$ and $s_{D,1}$, respectively. Indeed, a smooth transition is expected in the rigidity dependence of the diffusion coefficient, when the turbulence changes between two regimes [77, 78]. The second break at about 300 GV is directly observed in the AMS-02 data [22]. The fact that the break is more pronounced in secondaries than in primaries [79] clearly points to a change in the diffusion coefficient rather than in the injection spectrum [80]. A natural explanation for this scenario is for example provided by self-generated turbulence, as explored in Refs. [78, 77]. In contrast, a break at lower rigidities ($\sim 10$ GV) is more speculative because also other processes like convection and reacceleration influence the CR spectra [71, 23]. Therefore, we will explore two different models (as in Ref. [81]), one with the low-energy break and one without. The models will be detailed further below.

Convective winds with the velocity $\bm{V}_{\rm c}$ can transport the CRs away from the Galactic plane and induce adiabatic energy losses. We employ a constant convection velocity that is perpendicular to the Galactic plane, $\bm{V}_{\rm c}=v_{0,\rm c}\,{\rm sign}(z)\,\bm{e}_{z}$. Increasing the convection velocity decreases the CR flux below $\sim 10$ GV. In contrast, reacceleration is due to scattering of CRs on Alfvén magnetic waves and has the opposite effect on the CR spectra. In a head-on scattering with the Alfvén waves a CR gains energy while it loses energy in back-on collisions. However, head-on collisions are more likely such that statistically the CRs gain energy. In eq. (2), reacceleration is modeled as diffusion in momentum space with the coefficient $D_{pp}$. Larger magnetic turbulence makes reacceleration more efficient, such that $D_{pp}\sim 1/D_{xx}$ [82]. The amount of reacceleration further depends on the Alfvén velocity, $D_{pp}\sim v_{\rm A}$.

The term $\partial/\partial p(\mathrm{d}p/\mathrm{d}t\psi_{i})$ represents continuous energy losses from ionization and Coulomb collisions [83], while the last two terms of eq. (2) describe catastrophic losses by fragmentation and decay with the characteristic time scales $\tau_{f,i}$ and $\tau_{r,i}$, respectively.

Finally, when the CRs enter the heliosphere, they are deflected and decelerated by the solar winds. The effect is known as solar modulation and varies in a 22-year cycle. The effect is most prominent at low energies. In principle, solar modulation can be described by a diffusion equation, similar to equation (2), but with the geometry, magnetic field and turbulence adjusted to the heliosphere. There are semi-analytical [84] or fully numerical codes [85, 86] solving this equation numerically. We use a simplified approach instead, and treat solar modulation in the force-field approximation [87]. This approximates well enough the effect of solar winds for nuclei observed above about few GeV.

We explore two distinct frameworks for CR propagation labelled INJ.BRK and DIFF.BRK. The assumptions and free parameters of each model are detailed below.

The INJ.BRK (injection break) model has been explored extensively in literature [83, 88, 89, 90, 91, 92, 93]. It also matches the model that we studied in the previous work on antiproton constraints [27]. The assumptions and free parameters are sketched in figure 1 (upper row). We employ a broken power law for the injection spectra of the primary (astrophysical) CRs with slopes $\gamma_{1}$ and $\gamma_{2}$ below and above the break position at $R_{0}$ with a smoothing $s$. We use different slopes for $p$ and He in the fit to account for the observed difference of the slopes in the CR fluxes of the two primaries. Understanding this difference is subject to current research. Possible explanations are, for example, different source populations [94, 95, 96] or a $Z/A$-dependence of efficiency of Fermi-shock acceleration [97, 98]. The diffusion coefficient is modeled as a single broken power law with the break $R_{D,1}$ around 300 GV (i.e. $\delta_{l}=\delta$ in eq. (2.5)). Furthermore, we allow for reacceleration and convection⁴⁴4In the recent update of Galprop [99] from version 56 to 57 the implementation of convection was changed, correcting a false definition of the Crank–Nicolson coefficients near the Galactic plane. However, this has negligible impact on our results because our CR fits prefer very small (compatible with 0) convection velocities. with $v_{\rm A}$ and $v_{0,c}$ as the two free parameters, respectively. Finally, solar modulation is treated in the force-field approximation. We allow for a slightly different solar modulation potential for antiprotons because solar modulation is charge-sign dependent [100].

The DIFF.BRK (diffusion break) model employs a single power law for the injection spectrum of the primary CRs. This type of model has also been studied in literature (see e.g. [83]), but only more recently it is tested against AMS-02 data [101, 102, 23]. In contrast to the INJ.BRK model, reacceleration is replaced by an additional break, $R_{D,0}$, in the diffusion coefficient, see figure 1. It is conceivable that the interaction of CRs and magnetic turbulence causes a damping of turbulence at low energies that effectively leads to an increase of the diffusion coefficient [103].

For both setups, we fix the half-height of the diffusion halo to $4\,\mathrm{kpc}$, which is roughly the lower bound compatible with the beryllium data from AMS-02 [104, 105, 93, 23, 106]. This is different with respect to the propagation setup used in our previous work [27], in which $z_{h}$ was varied as free parameter. However, a precise determination of the halo size is currently prevented by large systematic uncertainty in the secondary fragmentation cross section [107, 104, 23]. Thus, at the moment, larger values for the half-height are equally viable because of the well-known degeneracy of $z_{h}$ with the normalization of the diffusion coefficient. We chose a small value for the half-height because this corresponds to the most conservative DM limit.

If we had performed the whole analysis with a different value of $z_{h}$, to a first approximation, we would have inferred a different value of $D_{0}$. The other propagation parameters would only change marginally, and also the astrophysical fluxes of primary and secondary CRs would not be affected by this. However, the DM flux would increase as a function of $z_{h}$. The difference between astrophysical CRs and CRs from DM annihilation is the spatial extent of the source term. The astrophysical CRs are produced in a thin layer of the Galactic plane, while DM annihilates and produces antiprotons in the entire diffusion halo [8]. Empirically, we found the following enhancement factor for the flux of DM antiprotons, which is valid for $z_{h}$ between 1 and 10 kpc:

	$\displaystyle f_{\bar{p},{\rm DM}}(z_{h})$	$\displaystyle=$	$\displaystyle 1.0+3.7\times 10^{-1}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm kpc}\right)+5.0\times 10^{-3}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm kpc}\right)^{2}$
			$\displaystyle-8.1\times 10^{-3}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm kpc}\right)^{3}+8.4\times 10^{-4}\left(\frac{z_{h}-4\,{\rm kpc}}{\rm kpc}\right)^{4}\,.$

We note that this factor is normalized to one at our benchmark of $z_{h}=4$ kpc.

The likelihood $\mathcal{L}_{\bar{p}}$ for the AMS-02 measurements of the antiproton flux can be written as $-2\log\mathcal{L}_{\bar{p}}=\chi^{2}_{\bar{p}}$, with

$$\chi^{2}_{\bar{p}}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})=\sum_{i,j}\left(\phi_{\bar{p},i}^{(\text{AMS})}-\phi_{\bar{p},i}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})\right)V^{-1}_{ij}\left(\phi_{\bar{p},j}^{(\text{AMS})}-\phi_{\bar{p},j}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})\right)\;,$$

(3.1)

where $\mathbf{\phi}_{\bar{p}}^{(\text{AMS})}$ denotes the AMS-02 measurements and $V$ is the covariance matrix. We estimate the covariance matrix following the approach of Refs. [71, 32], i.e. we write

$$V=\text{diag}(\sigma^{2}_{\text{stat},i})+\sum_{\alpha}V^{\alpha}+V_{\text{\rm xs}}\;.$$

(3.2)

The first two terms represent the uncertainty of the flux measurement by AMS-02, where $\sigma_{\text{stat},i}$ denotes the (uncorrelated) statistical uncertainty in bin $i$ and $V^{\alpha}$ denote additional systematic uncertainties.⁵⁵5AMS-02 provides only the diagonal entries of the sum of first two terms in eq. (3.2), i.e, it gives the uncorrelated statistical and systematic uncertainties in flux measurements. A crucial contribution to the correlated uncertainty comes from the antiproton absorption cross section in the AMS-02 detector, as carefully modeled in Ref. [32]. All other contributions can be approximated by

$$\left(V^{\alpha}\right)_{ij}=\Delta_{i}^{\alpha}\Delta_{j}^{\alpha}\exp\left(-\frac{(\log R_{i}/R_{j})^{2}}{2(l^{\alpha})^{2}}\right)\;$$

(3.3)

with $R_{i}$ denoting the rigidity of bin $i$, $\Delta_{i}^{\alpha}$ denoting the systematic error and $l^{\alpha}$ denoting the corresponding correlation length. Following Ref. [32] we consider errors in the effective acceptance, errors in the rigidity scale, unfolding errors, geomagnetic cut-off errors, template shape errors and selection errors.

In contrast to the first two terms, the third term takes into account the uncertainties in modeling the production of secondary antiprotons. Here we follow the procedure form Refs. [30, 31]. We translate the correlated uncertainties in the production cross section from Ref. [73] into the covariance matrix $V_{\rm xs}$ that can be applied directly on the antiproton flux. In more detail, we generate 2400 vectors of the cross section parameters using the correlated uncertainties from Ref. [73] and then calculate the antiproton flux with Galprop for each parameter vector. Finally, we have used the sample covariance for the antiproton flux evaluated at rigidities $R_{i}$ and $R_{j}$:

$$\left(V_{\rm xs}\right)_{ij}=\frac{1}{N-1}\sum\limits^{N}_{k=1}\big{(}\phi_{\bar{p},i}^{k}-\bar{\phi}_{\bar{p},i}\big{)}\big{(}\phi_{\bar{p},j}^{k}-\bar{\phi}_{\bar{p},j}\big{)}\;,$$

(3.4)

where the superscript $k$ denotes the parameter vector and $N=2400$ is the sample size. For this procedure, the CR propagation parameters are fixed to the best-fit values, meaning that $V_{\text{\rm xs}}$ is slightly different for the two propagation setups DIFF.BRK and INJ.BRK.

We want to gain an understanding of the reasonable parameter space of the parameters relevant for CR propagation. These will be our nuisance parameters in the likelihoods for the DM models we consider. We use MultiNest [108] as a tool to efficiently sample the parameter space and obtain the likelihoods of each sampled point based on our model. MultiNest uses a nested sampling algorithm that initially samples from the prior and then iteratively samples and stores parameter points from ordered prior volumes according to their respective likelihoods in the parameter space. We thus get an ensemble of samples covering the entire parameter space, with the most dense sampling within the regions of high likelihood.

For each propagation model described in section 2.2, we perform this fit twice. Our initial scan includes all of the propagation parameters described previously with the assumption of no DM signal contributing to the measured CR spectra, i.e. our null hypothesis. This fit results in our most general understanding of the parameter space. We choose the priors as uniform distributions that are either constrained by observations of different CR species or have been indicated by previous parameter inferences. We simultaneously fit the most recent AMS-02 antiproton-over-proton ratio together with the fluxes of protons and helium and the ${}^{3}\mathrm{He}/^{4}\mathrm{He}$ ratio [22, 109]. For proton and Helium we also use Voyager data [110]. For the MultiNest scans, we do not consider correlations in the data sets, i.e. we assume that the given error bars can be interpreted as uncorrelated statistical uncertainties. The correlations in the anti-proton data, modeled as discussed above, will be included when deriving DM constraints, see section 3.3.

In the fit, we consider the propagation parameters going into eq. (2) for each model. We consider the effect of solar modulation, which is modeled with the force-field potential $\varphi_{\mathrm{AMS-02}}$ for proton and helium. Antiprotons are allowed to have a different value of the force-field potential, as described earlier. In the fit this is taken into account by using the difference of the force-field potentials, $\varphi_{\mathrm{AMS-02}}-\varphi_{\bar{p}}$, as a free parameter. Furthermore, the normalization of the proton flux, $\mathrm{norm}_{p}$,⁶⁶6The global free normalization of all primary and secondary fluxes (except the DM flux) in Galprop is fixed by choosing the proton flux at 100 GeV. and the ⁴He isotopic abundance relative to the proton abundance, ${\rm Abd}_{{}^{4}\mathrm{He}}$, are free parameters. Finally, we allow for cross section uncertainties in the ³He production cross section by introducing nuisance parameters for the cross section normalization and slope, $A_{\mathrm{XS,^{4}\mathrm{He}\rightarrow^{3}\mathrm{He}}}$ and $\delta_{\mathrm{XS,^{4}\mathrm{He}\rightarrow^{3}\mathrm{He}}}$, as in Ref. [59]. The complete setup allows for independent constraints on the propagation and a fully Bayesian interpretation of all parameters.

The prior ranges and resulting best fit points of the scan (plus their 1$\sigma$ intervals) are shown in table 1. We show the result of the best fit parameter point applied to the fluxes of $p,\overline{p}$, ³He and He for both propagation models in figure. 2. While the best fit values for different physical properties vary significantly for both models, the overall fit to the data is very similar, in line with the findings of Ref. [59]. The combined $\chi^{2}$ here is 159.1 (163.7) for the DIFF.BRK (INJ.BRK) model. The total number of likelihood evaluations for the scan is $1.4\,(3.7)\cdot 10^{6}$.

In a second step, we also include a tentative DM signal and add the DM mass and annihilation cross-section as free parameters to the MultiNest scan. For simplicity, we fix the annihilation channel to $b\overline{b}$, which is a standard benchmark choice. Here, we are not primarily interested in the likelihood of the final result of the scan, but rather in how the posteriors of the propagation parameters shift in order to accommodate the additional DM signal. This is important for one of the two use-cases of the MultiNest samples for our results. In this work we consider more general DM models than the one considered in this scan, assuming annihilation only into $b\overline{b}$, but the overall shift in the posterior for the individual propagation parameters can well be approximated with this assumption. We have found in particular a shift for $\gamma_{1,p}$ and $\gamma_{1}$ at the level of about $1.7\sigma$ towards lower values, and for $\delta$ towards a higher value at about $1.5\sigma$ in the DIFF.BRK model. The remaining fit parameters, as well as the INJ.BRK model were affected less by the additional signal, i.e. the shift was at most $1\sigma$.

We use the combination of both scans to set up an extensive training set for the DarkRayNet simulation tool, a deep neural network emulator that we set up as in Ref. [27]. In the remainder of this work, this tool is used to speed up the required simulations of CR-spectra.

For both scans in each propagation model, we store all simulated spectra with a $\Delta\chi^{2}\leq$ 30 and combine them as a training set. We include model independent DM signals by sampling four sets of DM masses and annihilation branching fractions for each set of propagation parameters. We set a flat prior on a logarithmic scale for the DM mass between 5 GeV and 5 TeV, and consider annihilation into $q\overline{q},\ c\overline{c},\ b\overline{b},\ t\overline{t},\ W^{+}W^{-},\ ZZ,\ gg,\ hh$, as we expect non-negligible contributions to a resulting antiproton flux from these channels, and sample random branching fractions, always summing up to one. We don’t need to vary the annihilation cross section in the emulator, as it can be inserted as a normalization to the DM signal later. With the final training sample for each propagation model, the fully trained network can emulate accurate CR-spectra within an extended, relevant parameter space. With the addition of the second scan, the shift in preferred regions of the propagation parameters due to that additional signal are accounted for. More details on the DarkRayNet and the involved neural networks can be found in appendix A.1.

The second use-case of the MultiNest scans is the generation of a sample of propagation parameters, which can then be used to perform the marginalization as described in more detail in the following section. For this, we extract the equally weighted posterior samples from the scan.

The likelihood introduced above can be used to perform a simultaneous scan over DM and propagation parameters. In the context of constraining DM models, it is however often unnecessary to infer the preferred propagation parameters. In this case, it is convenient to eliminate the dependence on the propagation parameters and obtain a likelihood function that depends exclusively on the DM parameters. Since fixing the propagation parameters to specific values may lead to overly aggressive constraints, the two possible options are either profiling or marginalisation. The former approach, i.e. maximising the likelihood with respect to the propagation parameters for every choice of $\bm{x}_{\text{DM}}$ is not only computationally expensive, but may be sensitive to finely-tuned parameter regions where small excesses in the data can be fitted. Moreover, when using DarkRayNet this approach may be susceptible to inaccurate network predictions for parameter regions with insufficient training (see the more detailed discussion in Ref. [27]).

The most robust approach is therefore to calculate the marginalised likelihood

$$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})=\int\mathrm{d}\bm{\theta}_{\text{prop}}\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{\text{prop}})\pi(\bm{\theta}_{\text{prop}})\;,$$

(3.5)

where $\mathcal{L}_{\text{tot}}=\mathcal{L}_{\bar{p}}\cdot\mathcal{L}_{p}\cdot\mathcal{L}_{\text{He}}$ and $\pi(\bm{\theta}_{\text{prop}})$ denotes the assumed prior probability for the propagation parameters. In the following, we will adopt flat priors and consider the same parameter regions used for the MultiNest scans described above, i.e. we set $\pi(\bm{\theta}_{\text{prop}})=\Pi_{i}(\theta_{i,\text{max}}-\theta_{i,\text{min}})^{-1}$. In principle, the marginalised likelihood can be obtained through Monte Carlo integration by drawing a random sample $\{\bm{\theta}_{{\rm prop},i}\}$ of size $N$ from the prior probability distribution and writing

$$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})=\frac{1}{N}\sum_{\bm{\theta}_{{\rm prop},i}\sim\pi}\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm prop},i})\;.$$

(3.6)

In practice, a more accurate estimate⁷⁷7The gain in accuracy is both due to the fact that more weight is given to parameter regions where the likelihood is large, thus decreasing the uncertainty of the Monte Carlo integration, and due to the fact that DarkRayNet is only evaluated in parameter regions where sufficient training data is available to guarantee reliable predictions. is obtained by drawing samples from the posterior probability for the propagation parameters in the absence of a DM signal $\mathcal{P}(\bm{\theta}_{\text{prop}})\propto\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{\text{prop}})\pi(\bm{\theta}_{\text{prop}})$. In this case,

$$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})\propto\frac{1}{N}\sum_{\bm{\theta}_{{\rm prop},i}\sim\mathcal{P}}\frac{\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm prop},i})}{\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm prop},i})}\;.$$

(3.7)

The constant of proportionality is independent of $\bm{x}_{\text{DM}}$ and therefore drops out when calculating likelihood ratios.

Now we can make use of the fact that the contribution of DM to the local proton and Helium flux is completely negligible. Hence, the proton and Helium likelihoods $\mathcal{L}_{p,\text{He}}$ are independent of $\bm{x}_{\text{DM}}$ and therefore cancel in the likelihood ratio. We can therefore replace $\mathcal{L}_{\text{tot}}\to\mathcal{L}_{\bar{p}}$ in the equation above. The marginalised $\chi^{2}$ is then given by

	$\displaystyle\bar{\chi}^{2}(\bm{x}_{\text{DM}})$	$\displaystyle=-2\log\left(\sum_{\bm{\theta}_{{\rm prop},i}\sim\mathcal{P}}\frac{\mathcal{L}_{\bar{p}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm prop},i})}{\mathcal{L}_{\bar{p}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm prop},i})}\right)$
		$\displaystyle=-2\log\left[\sum_{\bm{\theta}_{{\rm prop},i}\sim\mathcal{P}}\exp\left(-\frac{\chi^{2}_{\bar{p}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm prop},i})-\chi^{2}_{\bar{p}}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm prop},i})}{2}\right)\right]$		(3.8)

with $\chi^{2}_{\bar{p}}$ defined in eq. (3.1). We note that by construction $\bar{\chi}^{2}(\bm{x}_{\text{DM}}=0)=0$.

In general, drawing a random sample from the posterior probability is far from easy. Fortunately, this sample only needs to be generated once and can then be used for arbitrary DM parameters. Moreover, such a sample is generated automatically as a by-product of the MultiNest scans that we have run to generate the training data for DarkRayNet and therefore requires no additional calculations. This sample of about $10^{4}$ sets of propagation parameters is included in pbarlike and used by default for marginalisation.

To conclude this discussion, we note that it is in fact possible to use a different definition of the antiproton likelihood for the marginalised likelihood than the one that has been used to calculate the posterior probability. In this case, eq. (3.8) simply becomes

$$\bar{\mathcal{L}}_{\text{tot}}(\bm{x}_{\text{DM}})\propto\frac{1}{N}\sum_{\bm{\theta}_{{\rm prop},i}\sim\mathcal{P}_{0}}\frac{\mathcal{L}_{\text{tot}}(\bm{x}_{\text{DM}},\bm{\theta}_{{\rm prop},i})}{\mathcal{L}_{\text{tot},0}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{{\rm prop},i})}$$

(3.9)

with $\mathcal{P}_{0}(\bm{\theta}_{\text{prop}})\propto\mathcal{L}_{\text{tot},0}(\bm{x}_{\text{DM}}=0,\bm{\theta}_{\text{prop}})\pi(\bm{\theta}_{\text{prop}})$. This means in particular that it is possible to modify the covariance matrix without the need to rerun the computationally expensive MultiNest scan.

Assuming all of DM annihilates only into bottom quarks, we calculate the marginalised log-likelihood ratios pbarlike, for each point on a 2D grid of DM model parameters $\{m_{DM},\langle\sigma v\rangle\}$. We consider the DM parameter ranges, $m_{DM}\in[10\ \mathrm{GeV},5\ \mathrm{TeV}]$ and $\langle\sigma v\rangle\in[10^{-28},10^{-24}]\ \mathrm{cm^{3}s^{-1}}$. An entire routine consisting of obtaining antiproton flux prediction from DarkRayNet, solar modulation and marginalised log-likelihood ratio calculation for a single point in DM model parameter space takes $\mathcal{O}(10)\ \mathrm{s}$. The evaluation of $\mathcal{O}(10^{4})$ antiproton fluxes (for the entire sample of $\mathcal{O}(10^{4})$ propagation parameter vectors used for marginalization) using DarkRayNet takes only $\mathcal{O}(1)\,\mathrm{s}$, whereas Galprop would take $\mathcal{O}(10)$   s for each set of propagation parameters. The most time-consuming part of the likelihood calculation turns out to be the modulation of each of $\mathcal{O}(10^{4})$ antiproton fluxes with appropriate solar modulation parameters.

The results are shown in figure 3. The region in colour correspond to $\bar{\chi}^{2}>0$ and hence prefers the case of antiprotons from a purely secondary origin. The region in gray corresponds to $\bar{\chi}^{2}<0$ and thus exhibits a small preference for the presence of a DM signal. The black line shows the $95\%$ CL upper bound for $\langle\sigma v\rangle$ as a function of DM mass $m_{DM}$ (corresponding to $\bar{\chi}^{2}=3.84$). For the DIFF.BRK propagation model (left panel), preference for DM signal is seen in $\mathrm{TeV}$ masses and for INJ.BRK (right panel) around $100\,\mathrm{GeV}$. The bound obtained using INJ.BRK model is similar to ones previously seen in literature [30, 31, 32, 24]. For the DIFF.BRK model, however, the DM preference is pushed to $\mathrm{TeV}$ masses. This can be understood with the help of residuals in figure 2. The additional break at low rigidities in the diffusion coefficient in the DIFF.BRK model allows for a good fit to low rigidity data, whereas it performs a poorer fit to data at high rigidities compared to INJ.BRK model.

To conclude this section, we note that pbarlike also returns log-likelihood ratios calculated using uncorrelated errors given by AMS-02. For more details, we refer to appendix B.

In the scalar singlet DM model, the SM is extended by a gauge-singlet real scalar boson $s$ stabilised by a $\mathbb{Z}_{2}$ symmetry [36, 37, 38]. The resulting scalar potential can then be written as

$$V(s^{2},H^{\dagger}H)=\lambda_{h}\left[\left(H^{\dagger}H\right)-\frac{v^{2}}{2}\right]^{2}+\frac{1}{2}\,\lambda_{hs}\,s^{2}\,H^{\dagger}H+\frac{1}{4}\,\lambda_{s}\,s^{4}+\frac{1}{2}\,m_{s,0}^{2}\,s^{2}$$

(4.1)

with $v=246\>\text{GeV}$ denoting the vacuum expectation value of the SM Higgs field. After electroweak symmetry breaking, we can replace $H=(h+v,0)/\sqrt{2}$ with a real scalar $h$ to obtain the potential

$$V(s^{2},h)=V(h)+\frac{1}{2}\,m_{s}^{2}\,s^{2}+\frac{1}{4}\,\lambda_{s}\,s^{4}+\frac{1}{2}\,\lambda_{hs}\,v\,h\,s^{2}+\frac{1}{4}\,\lambda_{hs}\,h^{2}\,s^{2}\;,$$

(4.2)

where $m_{s}=\left[m_{s,0}^{2}+\lambda_{hs}\,v^{2}/2\right]^{1/2}$ denotes the physical mass of the scalar singlet.

The quartic self-coupling $\lambda$ typically has no consequences for phenomenology, and we will therefore set it to zero in the following. The model is then fully characterised by the scalar singlet mass $m_{s}$ and the Higgs portal coupling $\lambda_{hs}$.⁸⁸8In the analysis below, we also include a number of nuisance parameters, most notably the local DM density $\rho_{\odot}$, which affects direct and indirect detection constraints. See Ref. [40] for details. The latter is responsible for all interactions of $s$ with SM particles, which determine the relic density of $s$ via the freeze-out mechanism and the potentially observable signatures of $s$ in satellites and laboratory experiments.

The phenomenology of the scalar singlet DM model has been discussed in great detail in Refs. [39, 40, 41]. Here, we repeat the analysis procedure from Ref. [40], i.e. we calculate the relic density $\Omega_{s}$ of scalar singlets and require $f_{s}\equiv\Omega_{s}/\Omega_{\text{DM}}\leq 1$, where $\Omega_{\text{DM}}=0.12$ [111]. Furthermore, we consider constraints from direct detection experiments (rescaled by the fractional abundance $f_{s}$ of scalar singlets), indirect detection experiments (rescaled by a factor $f_{s}^{2}$) and the LHC as well as the perturbativity requirement $\lambda_{s}<\sqrt{4\pi}$. In addition to previously considered constraints, we include the antiproton likelihood from AMS-02 discussed above, as well as a number of additional new likelihoods as detailed below.

In this section, we present updated constraints on the scalar singlet DM model. To do this, one could perform a comprehensive global scan with all the relevant available datasets and the resulting $\mathcal{O}(10)$ nuisance parameters. Instead, for a judicious use of computational resources, we take advantage of the publicly available results [112] (hereafter referred to as GC17) from the extensive global analysis performed by the GAMBIT collaboration in Ref. [40]. This analysis combined the then available relic density measurements, direct detection limits, limits from invisible Higgs decays, and indirect detection limits from dark annihilating into neutrinos and gamma-rays. The results in GC17 contain parameter spaces allowed by the limits considered and their corresponding combined likelihoods.

The large sample of $\mathcal{O}(10^{7})$ points in GC17 contains combined results from multiple sampling runs. These scans included the DM model parameters in the range $m_{s}\in[45\,\mathrm{GeV},10\,\mathrm{TeV}]$ and $\lambda_{hs}\in[10^{-4},10]$. For details on the nuisance parameters and their ranges, see table 2 of Ref. [40]. The scans identified the best-fit at $\lambda_{hs}=6.5\times 10^{-4},\ m_{s}=62.51\,\mathrm{GeV}$ corresponding to $f_{s}^{2}\langle\sigma v\rangle_{0}=8.65\times 10^{-28}\,\mathrm{cm^{3}s^{-1}}$.

For our analysis, we restrict the parameter space to $m_{s}\in[45\,\mathrm{GeV},5\,\mathrm{TeV}]$ to remain within the parameter region where DarkRayNet has been trained. We further reduce the GC17 sample to $\mathcal{O}(10^{5})$ points. We then update the old combined likelihoods in GC17 by adding the new likelihoods discussed above using the postprocessor scanner available within the GAMBIT sampling module ScannerBit [113]. We emphasize that this procedure is possible because there is no strong DM signal in the AMS-02 data, and hence we do not expect the new likelihoods to open up previously disfavoured parameter regions.

In a first step, we only add the new AMS-02 antiproton likelihood, to explore its impact on the allowed parameter space. The resulting updated constraints on the $f_{s}^{2}\langle\sigma v\rangle_{0}$–$m_{s}$ plane are shown in figure 4.⁹⁹9Plots were made using pippi [114]. Each panel shows the profiled likelihood in the DM model parameter space, normalized by the best-fit likelihood.¹⁰¹⁰10Note that the total likelihood here is profiled over the nuisance parameters in GC17 only; pbarlike already returns the marginalized likelihood and hence, the CR propagation and solar modulation parameters are not included in the GAMBIT postprocessing run as nuisance parameters. The contours show the $1\sigma$ and $2\sigma$ confidence regions, and the star indicates the best-fit point. The viable regions identified in previous studies, i.e. the resonance region with $m_{s}\approx m_{h}/2$ and small cross sections (left column) and the high-mass regions with $m_{s}\gg m_{h}$ and large cross sections (right column), are still retained. For comparison, see the right panel of figure 3 in Ref. [40].

The two rows of figure 4 correspond to the DIFF.BRK (top) and INJ.BRK (bottom) models. For both propagation models, the likelihood of the best-fit point is changed by $-2\Delta\log\mathcal{L}<0.1$. Neither propagation models shift the best-fit point to different masses. For DIFF.BRK, the best-fit point is however shifted to smaller couplings and a smaller log-likelihood, as this model leads to strong bounds on the rescaled annihilation cross section for DM masses in the 10–100$\,\mathrm{GeV}$ range (see also figure 3). For the INJ.BRK model, the small excess seen in figure 3 leads to a weaker bound on the rescaled cross section and a correspondingly larger coupling at the best-fit point with a larger log-likelihood. In the TeV range, on the other hand, the INJ.BRK model gives the stronger constraint than the DIFF.BRK model, leading to a smaller allowed parameter region. This observation highlights the relevant impact of the AMS-02 antiproton constraints on the scalar singlet DM model.

In a second step, we now also include the new direct detection constraints using DDCalc [5, 115] and the Higgs invisible branching ratio likelihood, which has been directly implemented within GAMBIT. The resulting updated constraints on the $\lambda_{hs}$–$m_{s}$ plane are shown in figure 5 for the DIFF.BRK (top) and INJ.BRK (bottom) models. As before, the figure also shows a zoom-in of the resonance region with mass in linear scale (left). The impact of the antiproton constraints is less visible in this case, as the relevant regions are already excluded by the new direct detection constraints. Hence, the top and bottom panels look almost identical. For both the propagation models, the likelihood of the best-fit is changed by $-2\Delta\log\mathcal{L}\approx 2.5$, driven by the small excess seen in Higgs invisible decays. Thus, the coupling of the best-fit is also shifted to larger values ($\sim 0.02$). In the resonance region, the allowed parameter space is constrained from the left by the improved direct detection constraints, which also exclude a large region of the parameter space in the high-mass region. The remaining parameter space will be further explored by future direct DM detection experiments [116].