Constraining the Higgs portal with antiprotons

The scalar Higgs portal is a compelling model of dark matter (DM) in which a renormalizable coupling with the Higgs boson provides the connection between the visible world and the dark sector. In this paper we investigate the constraint placed on the parameter space of this model by the antiproton data. Due to the fact that the antiproton-to-proton ratio has relative less systematic uncertainties than the antiproton absolute flux, we propose and explore the possibility to combine all the available $\bar{p}/p$ data. Following this approach, we are able to obtain stronger limits if compared with the existing literature. In particular, we show that most of the parameter space close to the Higgs resonance is ruled out by our analysis. Furthermore, by studying the reach of the future AMS-02 antiproton and antideuteron data, we argue that a DM mass of $\mathcal{O}(150)$ GeV offers a promising discovery potential. The method of combining all the antiproton-to-proton ratio data proposed in this paper is quite general, and can be straightforwardly applied to other models.


Introduction
Since the discovery of the first negative proton in 1955 [1], antiprotons have become a fundamental pillar in experimental high-energy physics, both in collider physics and astrophysics. In particular, antiprotons play a starring role in the context of indirect detection of dark matter (DM): they are copiously produced in the final stages of DM annihilations into Standard Model (SM) particles as a consequence of showering and hadronization processes. This is in particular true considering DM annihilation into hadronic channels, like for instance annihilation intobb. However, thanks to the electroweak radiative corrections, this is also true for DM annihilations into leptonic channels, providing that the DM mass is around the TeV scale [2][3][4].
Furthermore, the astrophysical background plaguing this potential DM signal is relatively well understood. In a standard scenario it mostly consists in secondary antiprotons originated from the interactions of primary cosmic-ray protons, produced in supernova remnants, with the interstellar gas. For these reasons the antiproton channel is considered one of the most promising probes to shed light on the true nature of DM [5][6][7][8].
However, all the experimental data collected so far show a fairly good agreement with the predictions of the astrophysical background, usually computed by means of dedicated codes such as GALPROP or DRAGON. Overturning the previous perspective, this negative results is often exploited to place strong bounds on the annihilation cross-section of DM into SM particles (see, e.g., refs. [9][10][11][12][13][14][15][16][17][18][19][20]). Following this line, in this paper we explore the constraining power of the antiproton data considering as a benchmark example the so-called Higgs portal DM model [21][22][23][24]. Despite its simplicity, in fact, this model offers a rich phenomenology, and it provides a simple and motivated paradigm of DM.
The aim of this paper is twofold. On the one hand, we extract our bounds focusing the analysis on the experimental data describing the antiproton-to-proton ratio instead of (as customary in the literature) the absolute antiproton flux. In this way we can get rid of the systematic uncertainties that usually preclude the comparison between data taken by different experiments. On the other one, we compare, in the context of the Higgs portal model and in a wide range of DM mass, our results with the bounds obtained considering the invisible Higgs decay width and the spin-independent DM-nucleon cross-section. The purpose of this comparison is to highlight the regions of the parameter space in which the antiproton data give the most stringent limits. This paper is organized as follows. In section 2, we briefly introduce the scalar Higgs portal model. In section 3, we discuss all the relevant aspects of our analysis; in particular, we present in detail the computation of thep p flux considering both the standard astrophysical background and the DM signal. In section 4 we present our results, and in section 5 we discuss future prospects. Finally, we conclude in section 6. In appendix A, we generalize our results to the fermionic Higgs portal model.

Setup: the scalar Higgs portal
The Lagrangian of the scalar Higgs portal model that we consider in this work has the following structure [21][22][23][24] where L SM is the SM Lagrangian, λ S is the Higgs portal coupling, and the real field S is a scalar gauge singlet with mass -after electroweak symmetry breaking -given by m S = (m 2 0 + λ S v 2 2) 1 2 ; H is the SM Higgs doublet with vacuum expectation value (vev) ⟨H⟩ = v √ 2 = 174 GeV.
The relevant parameter space of the model is the two-dimensional plane (m S , λ S ). After electroweak symmetry breaking the Higgs portal coupling generates the trilinear vertex L hS 2 = λ S vhS 2 2, where h is the physical Higgs boson; this interaction is responsible for all the phenomenological properties of the model since via this vertex the DM particle communicates with all the SM species. In this paper, we focus on the possibility that the scalar field S plays the role of cold DM in the Universe.
In appendix A, we will analyze a different type of Higgs portal with fermionic DM.

Relic density
Through the exchange of the Higgs in the s-channel, two DM particles can annihilate into all the SM final states that are kinematically allowed by the value of the DM mass, m S . The annihilation cross-section times the relative velocity v rel of the two annihilating DM particles takes the remarkably simple form where the square of the total energy in the c.o.m. frame is s = 4m 2 S (1 − v 2 rel 4). In eq. (2.2) Γ h ( √ s) is the off-shell decay width of the Higgs boson (with m * h = √ s), summed over all the SM final states. We use the public code HDACAY [25] to compute the width Γ h ( √ s). In this way we are able to include i ) the O(α s ) NLO QCD radiative corrections to the Higgs decay into quarks and ii ) the Higgs decay modes into off-shell gauge bosons. The importance of these radiative effects has been emphasized in ref. [26]. We plot in the left panel of fig. 1  s > 2m h = 252 GeV DM annihilation into two Higgses is kinematically allowed; the corresponding annihilation cross-section, however, cannot be recast in the form described by eq. (2.2) since, in addition to the s-channel exchange of the Higgs, also t-and u-channel diagrams in which the DM particle is exchanged contribute to the  Figure 1: Off-shell decay width of the Higgs boson Γ h ( √ s) in the energy interval √ s ∈ [40,200] GeV obtained using the public code HDECAY (left panel, solid lines), including the NLO QCD corrections to the Higgs decay into quarks and the Higgs decay modes into off-shell gauge bosons. We show separately the most important contributions relevant for the computation of the annihilation cross-section in eq. (2.2). The impact of the radiative corrections is proven by the comparison with the corresponding tree-level expressions (tree-level Higgs decay width into bb, dotted line, and cc, dot-dashed line); we also show the W + W − and ZZ kinematical thresholds (vertical dashed lines) to emphasize the importance of the Higgs decay modes into off-shell gauge bosons in the region close to the Higgs resonance, amplitude. We include this channel computing the cross-section analytically (see, e.g., ref. [27]). In eq. (2.2) Γ h,S represents the on-shell decay width of the Higgs boson and it consists of two pieces, namely Γ h,S ≡ Γ SM h + Γ h→SS ; Γ SM h = 4.217 MeV is the SM contribution while Γ h→SS is the decay width describing the process h → SS, kinematically allowed if m S < m h 2. The explicit expression of Γ h→SS is discussed in the context of the LHC bound (see section 2.2, eq. (2.8)).
The thermally averaged annihilation cross-section is given by where T is the temperature and K α=1,2 are the modified Bessel functions of second kind. We numerically solve the Boltzmann equation describing the evolution of the number density n(x) of the DM particles during the expansion of the Universe, being x ≡ m S T . In terms of the yield Y(x) = n(x) s(x), where s(x) is the entropy density, this equation reads GeV is the Planck mass, and g * (T ) is the number of relativistic degrees of freedom. At the equilibrium, we have where h eff (T ) is the effective entropy. 1 The integration of the Boltzmann equation gives the yield today, Y 0 , which is related to the DM relic density through The Planck collaboration has recently reported the value [30] Ω DM h 2 = 0.1199±0.0027 (68% C.L.); we compute the relic density according to eq. (2.7), and we impose to match the measured value. Before proceeding, let us stress that the value of the relative velocity sets the position of the Higgs resonance during DM annihilation. From eq. (2.2) it follows that the pole of the Higgs propagator is given by s = m 2 h ; at zero relative velocity, therefore, the resonant annihilation occurs at m S = m h 2 = 63 GeV, while in general it occurs at m 2 Considering the annihilation in the early Universe -i.e. taking for definiteness v rel = 1 2 -the resonance occurs at m S ≃ 61 GeV.

LHC bound
If m S < m h 2, the Higgs boson can decay into two DM particles resulting in the possibility to have a sizable invisible decay channel. In this case the invisible decay width of the Higgs is The invisible branching ratio is severely constrained by the current searches at the LHC [31, 32], and BR inv > 22% is excluded at 95% C.L. [33]. Using eq. (2.9) it is straightforward to translate this bound in the parameter space (m S , λ S ), and in section 4 we will include this constraint in our analysis. To give a quantitative idea of the size of the invisible branching ratio in eq. (2.9), notice that for the benchmark values m S = 20 GeV, λ S = 0.05 we have BR inv ≃ 73%. As soon as the invisible decay channel is kinematically allowed, therefore, it can easily dominate over the SM final states even for relatively small values of the Higgs portal coupling. Notice that invisible width in eq. (2.8) is equal to zero for the resonant value m S = m h 2. Therefore, it will be impossible to test this particular region using the bound on the invisible branching ratio.

Direct detection constraint
The Higgs portal interaction, through the exchange of the Higgs boson in the tchannel, provides the possibility to have a non-zero spin-independent elastic crosssection of DM on nuclei.
Integrating out the Higgs in the limit of negligible exchanged momentum, it is possible to write the following effective interactions between DM and light quarks and gluons inside the nucleus with q = u, d, s, and G 2 = G µν G µν , where G µν is the gluon field strength. Using this effective interaction, it is straightforward to compute the spin-independent crosssection describing the elastic DM-nucleon scattering where µ S ≡ m N m S (m N + m S ) is the DM-nucleon reduced mass, m N = 0.946 GeV is the nucleon mass, and f N = 0.303 is the Higgs-nucleon coupling [26]. The LUX experiment has reported the strongest limit on σ SI [34]. Using eq. (2.11) we translate this bound in the parameter space (m S , λ S ), and in section 4 we will include this constraint in our analysis. Before proceeding, let us stress a simple but important point. The square of the momentum transferred in a typical DM-nucleus elastic scattering always satisfies the condition −q 2 ≪ m 2 h , with q 2 = −2m Xe E rec where the mass of a Xenon nucleus is m Xe = 121 GeV and for the typical recoil energy one has E rec ∼ few keV. This simply implies that there is no resonant enhancement in elastic scatterings via the Higgs portal; as a consequence, the region with m S ≈ 63 GeV where the model reproduces the correct relic abundance is beyond the present reach of direct detection experiments, and it will be covered only in the next future.

Antiprotons: Background vs Signal
In this section we address the properties of the transport equation describing the propagation of cosmic rays in the Galaxy. In section 3.1, we discuss the background contribution, focusing our attention mostly on astrophysical background of protons and antiprotons. In section 3.2, we illustrate the strategy that we follow in order to extract our bound on the scalar Higgs portal model.

Selection of propagation models
Considering the total luminosity injected in the Milky Way galaxy via cosmic rays, most (∼ 90%) of it consists of primary protons, ∼ 10% of helium nuclei, a further ∼ 1% of heavier nuclei, and ∼ 1% of free electrons. In full generality, the evolution of the cosmic-ray density in the Galaxy is described by the following transport equation is the number density per total unit momentum of the ith atomic species, p is its momentum and β is its velocity. The construction of a propagation model consists in solving eq. (3.1) with a certain boundary condition for all the cosmic-ray species; in this way one can compute -for a given distribution and energy spectrum of the sources -the spatial distribution and energy spectrum after propagation. Eq. (3.1) contains a number of free parameters to be determined. Let us discuss these parameters one by one.
• For each nucleus, the source distribution and the injection index γ i (or indices, if a break is considered). In eq. (3.1) these informations are encoded into the source term Q i (p, r, z). Considering the contribution of the astrophysical background, Q i (p, r, z) describes the distribution and injection spectrum of supernova remnants [35]. As far as the spectral index is concerned, we assume -following ref.
[10] -the same spectral index γ i = γ for all the nuclear species.
• The normalization and energy dependence of the diffusion coefficient, D. In eq. (3.1), we assume the following functional form where ρ = pβ (Ze) is the rigidity of the nucleus of charge Z, D 0 is the absolute normalization at reference rigidity ρ 0 = 3 GV, δ is the diffusion spectral index related to the turbulence of the interstellar medium, and z t is the scale height that controls the vertical spatial dependence, which is assumed to be exponential; the halo thickness z h ≡ 2z t is the height of the propagation halo where stochastic diffusion and re-acceleration take place. An additional parameter, η, controls the low-energy behavior of the diffusion coefficient.
• The Alfvén velocity, v A . It parametrizes the efficiency of the stochastic reacceleration mechanism. In eq. (3.1) it enters in the explicit expression of the diffusion coefficient in momentum space, D pp [10].
• The convective velocity, ⃗ v c . It is the velocity of the convective wind, if present, that may contribute to the escape of cosmic rays from the Galactic plane. The convective velocity is zero in the Galactic plane and linearly increasing with the vertical distance z from the Galactic plane.
• For each nucleus, the scattering cross-sections on the interstellar medium gas.
The distribution of gas in the interstellar medium concentrates in the disk of the Galaxy. Its number density is denoted as n gas , which is mainly constituted by atomic and molecular hydrogen and helium. There are two main effects. On the one hand, the i-th nucleus is generated by the nuclear species j with cross-section σ ji ; on the other one, the i-th nucleus is destroyed by scattering on interstellar medium gas with total inelastic cross-section σ in i .
The propagation of cosmic rays can be simplified, if one takes into account only the high-energy region ( ≳ 10 GeV). In this regime, diffusion and energy losses play an important role while other effects, such as convection and re-acceleration, are negligible. However, following this approach one is forced to neglect all the cosmic ray data at low-energy; as a consequence, the ability of constraining DM modelsespecially with relative low mass -decreases dramatically. In order to take advantage of all the data set, the propagation equation needs to be solved numerically without these approximation. To achieve this result, we use the public code DRAGON [36,37], and our procedure goes as follows.
Following ref [10], we start from five different benchmark propagation models: KRA, KOL, CON, THK and THN. These propagation models are characterized by different halo height z t , slope of the diffusion coefficient δ, spectral index γ, and gradient of the convection velocity dv c dz. We collect the corresponding values in table 1. KRA, KOL and CON are characterized by the same halo height (z t = 4 kpc) but they describe differently turbulence effects and convection velocity. THK and THN, on the contrary, explore two extreme values for the halo height, namely z t = 10 kpc and z t = 0.5 kpc. Using these five benchmark propagation models, the intent is to capture a wide range of astrophysical uncertainties. For a more detailed discussion we refer the interest reader to ref. [10].
The second step in our analysis is to use the Boron-over-Carbon (B/C hereafter) data in order to determine -for each one of the propagation setups defined before -the remaining phenomenological parameters in eq. (3.1). The B/C data employed in our analysis come from the HEAO3 [38], CREAM [39] ATIC [40] and CRN [41] experiments. In table 1, we complete the definition of the five benchmark propagation models by minimizing the χ 2 against B/C data; in this way we obtain, as output of the fitting procedure, D 0 , η and the Alfvén velocity v A . The B/C data with energy larger than 0.5 GeV are considered, and solar modulation is fixed to the value Φ = 0.55 GV. The reason why we chose to define our benchmark propagation models focusing exclusively on the B/C data is that antiproton flux may originate both from astrophysical and exotic sources, while Boron and Carbon are usually generated only by astrophysical processes. Moreover, B/C represents the ratio between stable secondary cosmic-ray flux divided by the corresponding primary cosmic-ray flux, which is exactly the same asp p that we will use in the next step of our analysis. In the left panel of fig. 2 we show the best-fit value for the B/C flux for the five propagation models in table 1. Solid lines are obtained using Φ = 0 GV, and the dashed line is to tune the value of Φ in order to fit that low-energy data from the ACE collaboration. Notice that to compute χ 2 , we use fixed value of the solar modulation Φ = 0.55 GV to fit the other experiments data.
Finally, we can use the five propagation models in table 1 to compute the background contribution to thep p flux. The choice of thep p ratio is made with the purpose of decreasing the systematic uncertainties that come from the comparison of data taken by different experiments. As far as the antiproton flux is concerned, for instance, various experiments can have different absolute flux due to different energy calibrations, which we want to avoid. Using data describing thep p flux, on the contrary, we can safely combine different datasets. Thep p data employed in our analysis come from the BESS [42,43], CAPRICE [44] and PAMELA [45] experiments. At this stage, the only free parameter is the value of solar modulation; in principle, since different experiments operated at different time, we can use three independent values of Φ -one for each experiment -in order to fit the data, and in table 1 we show the result of the corresponding χ 2 fit (see caption for details).
Once the astrophysical background has been fixed, we are now in the position to discuss the contribution from DM annihilation in the Higgs portal model.

Antiproton bound on the scalar Higgs portal model
We now move to discuss the computation of the antiproton flux originated from DM annihilation in the context of the Higgs portal model.
We solve eq. (3.1) using as a source term where ρ DM is the DM density profile, dN dE p the antiproton emission spectrum, i.e. the number of antiprotons per each annihilation, and ⟨σv rel ⟩ 0 is the thermally averaged annihilation cross-section times relative velocity evaluated at the present epoch. The solution of eq. 3.1 allows us to compute -as a function of DM mass and portal coupling, and at the location of the Earth -the antiproton flux originated from DM annihilation, φ DM p (m S , λ S ). Next, we compute the localp p flux by combining DM  Table 1: Phenomenological parameters describing the five benchmark propagation models (first column) used in our analysis. In the next four columns of the table we collect the values of halo height z t , slope of the diffusion coefficient δ, spectral index γ, and gradient of the convective velocity dv c dz; these values are kept fixed, and define -for each propagation models -the corresponding properties of the diffusionloss equation (3.1). The normalization of the diffusion coefficient D 0 , the low-energy parameter η and the Alfvén velocity v A are obtained via a χ 2 fit of the B/C data from the HEAO3 [38], CREAM [39], ATIC [40] and CRN [41] experiments; in addition, we show the corresponding minimum χ 2 B C (divided by the number of degrees of freedom). Using the propagation models so defined, we compute the background contribution to thep p flux, and we fit the solar modulation potential Φ against data from the BESS, CAPRICE and PAMELA experiments. We show in the last two columns the corresponding values of χ 2 p p (against, respectively, the full dataset and the subset of PAMELA data). The reported values for the solar modulation potential Φ refer to the fit of the PAMELA data only. contribution and astrophysical background, Finally, by means of a χ 2 fit of thep p data, we extract a bound on the parameter space of the scalar Higgs portal model.
We will discuss in section 4.2 the impact of different DM density profiles on the results of our analysis. In the computation of the antiproton emission spectrum, we included -consistently with the computation of the relic density -the three-body final states consisting of one on-shell and one off-shell electroweak gauge bosons. Following ref. [46], we made use of the PYTHIA 8.1 event generator [47,48] to extract these energy spectra.

Results
Following the approach outlined in section 3, we derived the bound on the parameter space of the scalar Higgs portal model by analyzing the antiproton-to-proton ratio data, and in this section we present and discuss our main results. We show the bound as a 3-σ exclusion line in the planes (m S , λ S ) and (m S , ⟨σv rel ⟩ 0 ). In both cases we compare the antiproton bound with the region that reproduces the correct amount of relic abundance, according to the result of the numerical analysis outlined in section 2.1. In addition, we superimpose the constraints obtained considering the invisible Higgs decay width and the spin-independent DM-nucleon elastic crosssection as described, respectively, in sections 2.2 and 2.3.
In section 4.1 we analyze the impact of different propagation models while in section 4.2 we discuss the impact of different DM density profiles.

On the impact of different propagation models
In this section, we study the antiproton bound on varying the propagation model, according to table 1. We show our results in fig. 3 and in fig. 4 where, for definiteness, we use the NFW DM density profile [49] (see fig. 6, eq. (4.1), and table 2 below).
Let us start the discussion with some general comments. The region that reproduces the correct value of relic density is represented by a green strip, while the regions excluded by the LHC and LUX experiments are shaded, respectively, in purple and red. In fig. 3 we focus on small values for the DM mass, i.e. m S ∈ [25, 100] GeV, in order to emphasize the role of the antiproton bound in the region close to the Higgs resonance. In the left panel of fig. 3 we show our results in the parameter space (m S , λ S ); considering the green strip that reproduces the correct value of relic abundance, the resonant region is immediately recognizable because of the usual funnel-shaped form. As already discussed in section 2.1, this region extends also for values of the DM mass smaller than m h 2 = 63 GeV as a consequence of thermal effects during the freeze-out epoch, and this feature clearly emerges in the plot from the result of our numerical analysis; moreover, notice that both the bound on the invisible Higgs decay width and the spin-independent DM-nucleon cross-section can not rule The region above each one of these curves is excluded. The dashed gray line represents the bound placed by the Fermi-LAT experiment using the gamma-ray data from dwarf galaxies (see text). Left panel. Bounds on the parameter space (m S , λ S ). Right panel.
Bounds on the thermally averaged annihilation cross-section times relative velocity at zero temperature.
out this region because of the kinematical reasons discussed in sections 2.2, 2.3. In the right panel of fig. 3, we translate our analysis in the plane (m S , ⟨σv rel ⟩ 0 ). Away from the resonance the value of ⟨σv rel ⟩ 0 that reproduces the observed relic abundance is close to the usual WIMP-miracle cross-section ⟨σv rel ⟩ 0 ≈ 2 × 10 −26 cm 3 s −1 .
Close to the resonance, on the contrary, this value is distorted by the presence of the aforementioned thermal effects that move the position of the resonance during the freeze-out epoch. In particular, we find that for 50 ≲ m S ≲ m h 2 the thermally averaged annihilation cross-section times relative velocity today can be as small as ⟨σv rel ⟩ 0 ≈ 10 −29 cm 3 s −1 , while for m S ≈ m h 2 we have ⟨σv rel ⟩ 0 ≈ 10 −22 cm 3 s −1 . In fig. 4 we focus on large values for the DM mass, i.e. m S ∈ [100, 3000] GeV. As in fig. 3, we show the constraints placed by our phenomenological analysis in the plane (m S , λ S ), left panel, and (m S , ⟨σv rel ⟩ 0 ), right panel. In fig. 3 and fig. 4 we show the antiproton bound that corresponds to the five propagation models defined in section 3, and we use the same color code introduced in table 1. The results of our analysis point towards the following remarks.
• On a general ground, we find that the KOL, THK, CON and KRA propagation models give similar bound, while the THN propagation model places the weaker constraint. In greater detail, in the low mass region the KOL model, characterized by a strong re-acceleration, gives the strongest constraint. On the contrary in the high mass region the CON model, characterized by a strong convective wind, provides the most stringent bound. The presence of strong convective effects, in fact, hardens the antiproton flux thus leading to stronger constraints on heavier DM models [10]. It is worth noticing that the THN model, based on a thin diffusion zone, is disfavored by recent studies on synchrotron emission, radio maps and low energy positron spectrum [50].
• In the low mass region the antiproton bound is competitive with the bound obtained from direct detection and invisible Higgs decay width. In particular, we find that the bound from antiproton is the only one able to rule out the resonant region with m S ≈ m h 2. Let us stress once again that this specific value for the DM mass would be otherwise inaccessible. On the one hand, in fact, the invisible Higgs branching ratio goes to zero moving towards the kinematical threshold m S = m h 2; on the other one, the square of the momentum transferred in a typical DM-nucleus elastic scattering always satisfies the condition −q 2 ≪ m 2 h , with q 2 = −2m Xe E rec where the mass of a Xenon nucleus is m Xe = 121 GeV and for the typical recoil energy one has E rec ∼ few keV.
• In the high mass region the antiproton bound obtained using the KOL, THK, CON and KRA propagation models is competitive with the exclusion curve traced by the LUX experiment. In particular, as clear from the right panel of fig. 4, using the KOL, THK and CON propagation models it is possible to probe the thermal cross-section up to m S ≈ 160 GeV.
• For comparison, we show in the right panel of figs. 3, 4 the 95% C.L. exclusion curve obtained considering the measurement of the gamma-ray flux from the dwarf spheroidal satellite galaxies of the Milky Way [51]. These dwarf galaxies are some of the most DM-dominated objects known, and -because of their proximity, high DM content, and lack of astrophysical backgrounds -they are usually considered to be the most promising targets for the indirect detection of DM via gamma rays. For simplicity, in our analysis we used only the data from the Draco dwarf spheroidal galaxy since it gives the strongest constraint.
Let us now describe in more detail our approach. In order to use the result of ref. [51], first we computed the gamma-ray flux from the Draco dwarf spheroidal galaxy in the Higgs portal model under scrutiny, combining all the different annihilation channels including three-body final states. Then, for each value of the DM mass, we compared the gamma-ray flux previously obtained with the 95% C.L. exclusion limit in each of the 24 energy bins analyzed in fig. 2 of ref. [51]. Finally, we extracted the bound on the cross-section from the energy bin that provides the strongest constraint. We find that, both in the low-and high-mass regions, the bound from antiproton that we obtain using the KOL, THK, CON and KRA propagation models is more than one order of magnitude stronger than the bound obtained from the analysis of the gammaray flux measured from the Draco dwarf spheroidal galaxy. Needless to say, a more detailed analysis would require to include all the 25 dwarf spheroidal galaxies studied in ref. [51] together with a more careful investigation of the systematic errors involved. This task goes well beyond the purpose of the simple estimation that we derived in this work, and will be left for future investigation.
In conclusion, we have found that the antiproton bound provides a strong constraint on the parameter space of the scalar Higgs portal model introduced in section 2. Remarkably, the constraining power of the antiproton data is comparable to the exclusion curves placed by the LHC and LUX experiments in particular for m S ≳ 50 GeV. Most importantly, the antiproton bound is the only one able to rule out the Higgs resonant region for m S ≈ 63 GeV. This conclusion does not strongly depend on the model used to describe the dynamics underlying the propagation of charged particle in the Galaxy; in particular, we have shown that the KOL, THK, CON and KRA propagation setups give, in magnitude, similar bounds. In order to stress this point we show in fig. 5, following ref. [52], a zoom on the Higgs resonant region. We introduce the variable ∆ ≡ (2m S − m h ) m h , and we present our constraints in the plane (∆, λ S ). From this point of view, it is clear that the bound we get from the antiproton data is by far the most stringent if compared with LHC and LUX results. The only region left unconstrained by the antiproton bound is the small mass window with 10 ≲ −1 ∆ ≲ 10 3 which corresponds to 56.7 ≲ m S ≲ 62.9 GeV with 10 −4 ≲ λ S ≲ 10 −2 . In this region the position of the Higgs resonance is subject to thermal effects; as previously discussed, in this small mass window the resonant annihilation cross-section reproducing in the early Universe the correct relic abundance corresponds to an off-resonant value in today's annihilations in the Galactic halo.
In this section we extracted the antiproton bound using the standard NFW profile in order to describe the density distribution of DM in the Galaxy. In the next section, we will discuss the impact of different DM halo profiles.

On the impact of different DM density profiles
In this section we explore the impact of different DM density profile on the analysis of the antiproton data. In addition to the NFW profile [49] already used in section 4.1, we repeat our analysis using the Einasto [53,54] and the Isothermal profile [55]. The former -similar to the NFW profile and characterized by a DM density distribution peaked towards the Galactic center -is favored by the latest standard numerical simulations [56,57] while the latter -characterized by a constant core -seems to be in agreement with the numerical simulations that include baryons, because of large exchange of angular momentum between the gas and DM particles [58]. We show these three DM density distributions in fig. 6 and eqs. (4.1, 4.2, 4.3), while in table 2 we collect the numerical values of the parameters that enter in their definitions.   fig. 6. For the Einasto profile, α = 0.17.
We show our results in fig. 7, for the low mass region, and in fig. 8, for the high mass region. In both cases we focus on the plane (m S , ⟨σv rel ⟩).
Let us start our discussion pointing out an important argument to keep in mind for the rest of the section. As a rule of thumb, one would expect that the antiproton flux from DM annihilation is larger for profiles in which the DM density is enhanced towards the Galactic center while is smaller for density distribution described by an isothermal sphere; as a consequence, one would naïvely guess that a common feature of the analysis is that the antiproton bound is always more (less) stringent for the NFW (Isothermal) profile. In general, however, this conclusion turns out to be partially incorrect. What really matters, in fact, is not the value of DM density at the Galactic center but at the position where the antiprotons -whose flux is measured on Earth -are generated. As one can easily imagine, further insights on this issue are strongly linked to the assumed propagation model, and our analysis points towards the following results.
• We find that the antiproton bound obtained using the THN model, based on a very thin diffusion zone, does not significantly depend on the DM density profile, and we do not show the corresponding plot in figs. 7, 8. In the THN model, in fact, the antiproton flux from DM annihilation is dominated by local contributions where the three profiles are equivalent [10].
• As far as the antiproton bound obtained assuming the CON propagation model is concerned, we find that also in this case the impact of different DM density distribution is moderately negligible, as shown in the upper-left panel of figs. 7, 8. As already observed in ref. [10], therefore, we argue that the uncertainty related to the DM distribution towards the Galactic center has a negligible effect in the CON model in which the antiproton flux from DM an- nihilation is dominated by local contributions.
• The impact of different DM density distribution is relevant considering the antiproton bound obtained using the KOL, THK and KRA propagation models. In these models, therefore, a large contribution on the antiproton flux from DM annihilation comes from non-local regions pointing towards the Galactic center where the three profiles present sizable differences, being more or less peaked. In greater detail, we find that the Isothermal (Einasto) profile gives the weaker (stronger) constraint; moreover, comparing with the NFW case, the bound obtained assuming the Isothermal density distribution has the largest deviation, while the Einasto density distribution gives a similar result. Comparing the three profiles, as done in fig. 6, we notice that in the region r ≳ 0.5 kpc -being r the radial distance from the Galactic center -the density distribution in both the Einasto and the NFW profiles are significantly larger than the Isothermal case; moreover, in this region the Einasto density distribution is larger if compared with the NFW one, thus reflecting the hierarchy observed in the exclusion curves. On the contrary, for r ≲ 0.5 kpc, the density distribution in the NFW case is larger w.r.t. the Einasto profile. All in all, we argue that in the KOL, THK and KRA propagation model the antiproton flux from DM annihilation is dominated by regions close to the Galactic center, with 0.5 ≲ r ≲ r ⊙ kpc. In order to strengthen this argument, we show in fig. 9 the local antiproton flux coming from DM annihilation and, in the inset plot, the relative contribution from a region enclosed within 1 kpc from the Galactic center. For definiteness, we consider ⟨σv rel ⟩ 0 = 3×10 −26 cm 3 s −1 and DM annihilation intobb (W + W − ) for m S = 70 GeV (m S = 700 GeV). 2 From this plot it is clear that for all the propagation models the contribution to the total antiproton flux from the inner Galactic region is at most 20%. Moreover, as expected, the THN and CON propagation models receive a negligible contribution from the region with r < 1 kpc. The thin height of the diffusion zone and the presence of strong convective wind efficiently remove a large fraction of the antiprotons originated towards the center of the Galaxy increasing their escape probability.
In conclusion, we have shown that the role of the DM density distribution in the Galaxy plays only a relatively marginal role in our analysis, and the astrophysical uncertainties affecting the limits on the scalar Higgs portal model that one can obtain using the antiproton data are mostly dominated by the details of the propagation model. 3 In this regards, the antiproton data that will be released by the AMS-02 experiment will play a crucial role in order to improve the current sensitivity. In the next section, therefore, we will briefly discuss future perspectives.  Figure 9: Local antiproton flux from DM annihilation for two representative values of DM mass and the five propagation models defined in table 1. In the inset plots, we show the relative contribution to the total local flux coming from a region enclosed within 1 kpc from the Galactic center. We use the NFW profile, and we do not include solar modulation effects, Φ = 0 GV.

Future perspectives
In this section we discuss some future perspectives related to our analysis. In section 5.1, we analyze the constraining power of the antiproton data that will be released by the AMS-02 experiment. In section 5.2, we analyze the antideuteron flux from DM annihilation in the scalar Higgs portal model.

AMS-02
The Alpha Magnetic Spectrometer (AMS-02) is a particle physics detector hosted on board of the International Space Shuttle, and designed to measure various cosmic-ray fluxes; thanks to this instrument, in the next future a more precise determination of the antiproton and antiproton-to-proton ratio will improve the constraints derived in this paper. To get a more concrete idea, we can simulate the prospects of this experiment by means of a set of mock data, in the energy range of (1 GeV, 450 GeV). Following ref. [59], we use a linear approximation of the AMS-02 detector energy resolution ∆E E = 0.042 E GeV + 10 % , (5.1) [11,18]. which determines the energy bin-size of the data. Having the bin size ∆E i and the flux Φ i for each bin, the observed antiproton number can be derived as Np = ap Φ i ∆E i ∆t, where we take ≃ 1 for the efficiency, ap = 0.2 m 2 sr for the geometrical acceptance of the instrument, and ∆t = 1 year for the reference data taking time. Since the dominant statistical uncertainty comes from the antiproton rather than the proton flux, the statistical error is approximately 1 Np; for definiteness, we fixed the systematic uncertainties to be 5% for one-year data taking. The uncertainty at each data point is the sum in quadrature of the systematic and statistical errors. The central value of the data for each bin, Φ i , follows the predicted flux from our five benchmark propagation models, which does not contain any DM contribution. With these mock data in hand, we can study the future sensitivity of antiproton-to-proton ratio data on DM models. In fig. 10 we show the projected bounds on the DM annihilation cross-section in the low-(left panel) and high-mass region (right panel) considering, as two extreme cases, the THN and KOL propagation models. We find that the future AMS-02 antiproton data will improve the bound for more than an order of magnitude in the high-mass region. On the contrary, in the low-mass region, our analysis show a little improvement if compared with the existing data. This is because we already exploited the full available data set, which is of reasonably good quality in the lowenergy region, and goes until the lowest value of 0.1 GeV, while the mock AMS-02 data starts from 1 GeV. In the high-energy region, on the contrary, AMS-02's resolution and luminosity are much better than the current status.

Antideuteron
Antideuteron has been proposed in ref. [60] as a promising indirect signal of DM annihilation in the Galactic halo. On a general ground, the annihilation of DM particles into SM hadronic channels -i.e. qq, W + W − , and ZZ -may produce antideuterons in the final state as a consequence of a two-step process. First -after showering, hadronization, and decay of unstable particles -a large number of antiprotons and antineutrons are produced. Second, an antiproton-antineutron pair may coalesce to form an antideuteron nucleus. The description of this process is usually addressed in the context of the so-called coalescence model. Given an antiproton and an antineutron with four-momenta k μ p and k μ n , the coalesce model approximates the probability for the formation of an antideuteron with the step function Θ(∆ 2 + p 2 0 ), where ∆ µ ≡ k μ p −k μ n , and p 0 is the maximum value of relative momentum that allows to form an antideuteron bound state. In this picture, p 0 is a free parameter and its numerical value has to be extracted from experimental data (see refs. [61,62] for a detailed discussion). In our analysis we use the results of ref. [63] in order to reconstruct the antideuteron energy spectra produced by DM annihilation. These energy spectra have been obtained using p 0 = 160 MeV, and the coalescence model previously discussed has been applied studying DM annihilations event-by-event [64]. Let us now move to discuss the antideuteron produced by high-energy astrophysical phenomena. The most relevant argument supporting the claim in ref. [60], in fact, emerges from the comparison between the antideuteron signal produced by DM annihilation and the corresponding astrophysical background. To be more precise, there are two key points to keep in mind. i ) Antideuterons are produced by high-energy collisions between extragalactic cosmic rays (mostly p,p and He) and the interstellar gas (mostly H and He) in the Galactic disk; the corresponding cross-section is small and -most importantly -it is characterized by a relatively high kinematical threshold; for instance the energy threshold for the creation of an antideuteron from a collision of a cosmic ray proton (antiproton) with the interstellar gas is E th = 17 m p (E th = 7 m p ). Let us give a closer look to this numbers considering for definiteness the scattering between a cosmic ray proton and a proton at rest in the interstellar gas. Because of conservation of baryon number, the production of an antideuteron from a protonproton collision requires a six-body final state, with a total energy squareŝ > (6m p ) 2 . On the other hand, in the rest frame of the gas,ŝ = (m p + E p ) 2 − k 2 p , where E p and k p are the energy and momentum of the cosmic ray proton. From these considerations, it follows that the impinging cosmic ray proton needs to have an energy E p > 17 m p in order to create an antideuteron. ii ) The binding energy for an antideuteron is extremely low, namely Bd ≈ 2.2 MeV. This implies that antideuterons are easily destroyed, and -as a consequence -they do not have to possibility to propagate long enough in order to loose most of their energy. The astrophysical background of antideuterons with kinetic energy Ed ≲ 1 − 3 GeV, therefore, is expected to be extremely small. Below these kinetic energies, an antideuteron flux originated from DM annihilation may easily stand out from the astrophysical background for more than one order of magnitude. In order to translate these qualitative statements into more quantitative results, we need to compare antideuteron background and DM signal after propagation in the Galaxy. The cross-sections describing production, elastic and inelastic scattering of antideuteron -key ingredients in order to solve the corresponding propagation equation, see eq. (3.1) -are not well known. Following ref. [65], we implemented in the DRAGON code all these cross-sections using the results of ref. [66], where they were extrapolated from experimental data under reasonable assumptions. In this way, we are in the position to solve numerically the propagation equation for antideuterons considering both the astrophysical background and the DM signal.
Present and future experiments -two energy bands of the AMS-02 experiment and the three phases of the General AntiParticle Spectrometer (GAPS) [67,68] -will increase the sensitivity of searches for cosmic-ray antideuteron over the current bound set by the BESS experiment. For the proposed sensitivities of AMS-02 and GAPS experiments we use the values from ref. [65]. In fig. 11, we present the predictions of the scalar Higgs portal model for the antideuteron flux, comparing background and background plus DM signal hypothesis. We analyze two benchmark values for the DM mass, and we use for definiteness the NFW density profile. In the left panel of fig. 1 we take m S = 160 GeV, with ⟨σv rel ⟩ 0 = 2×10 −26 cm 3 s −1 . These values correspond to a DM candidate that, according to figs. 4, 8, lies close to the present bound placed by the analysis of the PAMELA antiproton data. We find that the corresponding total antideuteron flux is higher than all the experimental sensitivities of the GAPS experiment assuming the KOL, THK and KRA propagation models. Therefore in these cases, if such DM candidate is realized in Nature, we expect -in principlea combined detection in both antideuteron and antiproton channels. For illustrative purposes, we show in the right panel of fig. 11 a different situation with m S = 400 GeV, ⟨σv rel ⟩ 0 = 2 × 10 −26 cm 3 s −1 . These values correspond to a DM candidate that, according to fig. 10, lies close to the future bound that will be placed by the analysis of the AMS-02 antiproton data. In this case, however, the total antideuteron flux will be hardly distinguishable from the astrophysical background.

Conclusions and outlook
In this work, we extracted a new bound on the scalar Higgs portal DM model using high-energy cosmic-ray astrophysics. In summary, the main points of our analysis are the following. First, we studied the propagation equation that governs the motion of charged particles in the Milky Way galaxy using the numerical code DRAGON. This equation depends on several free parameters that encode the astrophysical uncertainties describing the propagation of cosmic rays in the Galaxy. We fixed some of these parameters -i.e. the halo thickness, the source spectral index, the rigidity slope and the gradient of the convection velocity in the vertical direction -so to define five different propagation setups. Using the measurement of the boron-to-carbon ratio performed by the HEAO-3, ACE, CREAM, ATIC and CRN experiments, we fixed, via a minimization procedure, the remaining ones -i.e. the normalization of the diffusion coefficient, the Alfvén velocity and the low-energy diffusion index. Second, using the propagation setups fixed by the B/C data, we predicted the background contribution to the antiproton-to-proton ratio. Finally, we computed the antiproton flux from DM annihilation in the Higgs portal model including three-body final states and QCD radiative corrections; using the antiproton-to-proton ratio previously discussed, and combining background and DM signal, we extracted 3-σ bound on the parameter space of the model. The use of the antiproton-to-proton ratio allowed us to combine different experiments, namely PAMELA, BESS and CAPRICE data. In the antiproton-to-proton ratio, in fact, several systematic effects that plague the comparison between different experiments -as for instance different energy calibration -are integrated out. We compared our antiproton bound with the constraints coming from the LUX and LHC experiments considering -respectively -direct de-tection of DM particles and the invisible decay width of the Higgs boson. At the same time, we required to reproduce the observed amount of relic density.
We found that the antiproton bound is competitive; in particular, it provides the most stringent constraint on the model in the mass range m S ≈ 80−300 GeV for most of the analyzed propagation setups. Most importantly, the antiproton bound is the only one able to put in significant tension the resonant region m S ≈ m h 2, otherwise of difficult access to direct detection and collider searches.
In our analysis, we investigated the impact of astrophysical uncertainties related to different propagation setups and different models for the DM density distribution in the Galaxy. Moreover, we discussed future perspectives using a set of mock data in order to simulate those that will be released in the near future by the AMS-02 experiment. Finally, we highlighted in the context of the scalar Higgs portal model the role of the antideuteron channel as an important indirect detection observable able to provide a signature of annihilating DM.