Dark matter in the Higgs resonance region

The singlet scalar Higgs portal model provides one of the simplest explanations of dark matter in our Universe. Its Higgs resonant region, $m_\text{DM}\approx m_h/2$, has gained particular attention, being able to reconcile the tension between the relic density measurement and direct detection constraints. Interestingly, this region is also preferred as an explanation of the Fermi-LAT $\gamma$-ray Galactic center excess. We perform a detailed study of this model using $\gamma$-ray data from the Galactic center and from dwarf spheroidal galaxies and combine them with cosmic-ray antiproton data from the AMS-02 experiment that shows a compatible excess. In the calculation of the relic density, we take into account effects of early kinetic decoupling relevant for resonant annihilation. The model provides excellent fits to the astrophysical data either in the case the dark matter candidate constitutes all or a subdominant fraction of the observed relic density. We show projections for future direct detection and collider experiments to probe these scenarios.


I. INTRODUCTION
The origin of dark matter (DM) is one of the main unsolved puzzles in fundamental physics today, implying physics Beyond the Standard Model (BSM).The Singlet Scalar Higgs portal (SHP) model is among the simplest UV complete BSM theories that provide a valid DM candidate [1][2][3][4][5][6][7][8][9].While a large portion of the model parameter space has already been excluded by the interplay of direct detection and relic density constraints, the Higgs resonant region, m DM ≈ m h /2, is particularly challenging to probe.It constitutes one of the two remaining viable regions [10,11] within the model.Independently, a DM mass of around 60 GeV is also preferred as an explanation of two potential indirect detection signals of DM [12] as we detail below.In fact, putting together all relevant observational data allows for a scenario where the DM candidate constitutes all or just a subdominant fraction of the observed DM density and reveals an intricate relation between the model parameters and the DM fraction.Due to the sharp resonance, the results are highly sensitive to the DM mass in the region m DM = m h /2±O(Γ h ), where Γ h denotes the total Higgs width.Therefore, the Higgs resonant region deserves a closer look, which we will provide in this paper.
Several groups have reported the detection of an excess of γ rays, which is labeled as the Galactic center excess (GCE), in the data of the Fermi Large Area Telescope (Fermi-LAT) in the direction of the center of the Milky Way (see, e.g., [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]).Recently, Refs.[27,28] have provided comprehensive and updated results for the GCE properties.They find that the GCE Spectrum Energy Distribution (SED) is peaked at a few GeV and has a high energy tail significant up to about 50 GeV.The spatial distribution of the GCE is compatible with a DM template modeled with a generalized Navarro-Frenk-White density profile with slope γ = 1.2−1.3.Moreover, the GCE centroid is located at the dynamical center of the Milky Way and its morphology is roughly spherically symmetric and does not change with energy.Therefore, all the characteristics of the GCE are perfectly compatible with γ rays produced from DM particles annihilating in the main halo of the Milky Way.The GCE SED can be well modeled with DM particles of mass (40−80) GeV annihilating into b b with a thermal annihilation crosssection [12,20,21,24,29], which is the correct crosssection to explain the observed relic density of DM in the Universe [30].Among several BSM theories that have been proposed to fit the GCE, the SHP model with a DM mass close to the Higgs resonance provides one of the simplest solutions (see, e.g., [12,24,[31][32][33][34]).
An alternative interpretation is associated to a population of millisecond pulsars located around the Galactic bulge that would produce a signal with properties compatible with the GCE (see, e.g., [35][36][37][38]).Outbursts of cosmic rays (CRs) from the Galactic center have also been proposed as possible interpretations for the GCE (see, e.g., [39][40][41]).
In addition to γ rays, other cosmic particles, i.e. messengers, could be produced from DM particle annihilation, such as antiprotons (p).Interestingly, different groups have found an excess of p, with respect to the secondary production, in the data of AMS-02 [53] between (5 − 20) GeV.Its significance has been found to vary from 1 to 5σ depending on the analysis technique, CR propagation model and, more importantly, whether an estimation of the correlations in the AMS-02 systematic errors is considered [12,29,[54][55][56][57][58][59] (see, e.g., [60] for a review).DM particles with a mass of (60−80) GeV annihilating into b b can explain this excess, and a possible link to the GCE has been investigated [12,54,56,57].
Laboratories on Earth are also sensitive to BSM particles which could form DM. Experiments located at LHC, the largest operating particle collider to date, test BSM physics, including DM, making use of a variety of production processes and final state signals (see, e.g., [61] for a review).The LHC run 3 that started in July 2022 will further increase the sensitivity to DM signals.
Direct detection experiments such as LUX-ZEPLIN [62] (LZ) and XENONnT [63] are improving significantly our discovery potential of interactions between DM particles and detector atoms.Both experiments have recently published the first results for the upper limits of the spin-independent cross sections for 60 days and 97.1 days of data taking, respectively, setting new limits by about a factor of two stronger with respect to previously existing measurements.These experiments will push the limits closer to the neutrino fog [64].Next generation direct detection experiments such as DARWIN [65] will be able to investigate cross sections a factor of about 10 smaller than LZ and XENONnT (see, e.g., [66] for a review).
In this paper, we investigate the DM interpretation of the GCE within the SHP.We use a multi-messenger and multi-strategy approach.In particular, we perform a combined analysis of different cosmic messengers such as γ rays, detected from dSphs and the Galactic center, and the flux data of p from AMS-02.Then, we take advantage of the complementarity between direct and indirect detection and collider searches to verify whether the DM properties that explain the cosmic flux data are compatible with the LZ and LHC constraints.Finally, to find the model parameter space that matches all the observations, we compute the relic density, taking into account effects of kinetic decoupling that are relevant for resonant DM annihilation [67,68] during freeze-out.
The paper is organized as follows.In Sec.II, we intro-duce the SHP model and compute the annihilation cross sections and spectra relevant for indirect detection.In Sec.III, we calculate the DM relic density.Constraints from collider and direct detection experiments are discussed in Secs.IV and V, respectively, before interpreting γ-ray and p data within the model in Sec.VI.Finally, in Sec.VII, we combining direct, indirect, collider and cosmology data discussing implications for the viable parameter space.We conclude in Sec.VIII.A Python package, called SingletScalar DM, enabling the fast calculation of the dark matter spectra, relic density, as well as direct detection and collider constraints within the model is available on Github.

II. SINGLET SCALAR DARK MATTER WITH
A HIGGS PORTAL

A. Overview of the model
The SHP model is among the simplest DM models extending the SM by just a DM particle candidate, S, taken to be a real scalar and a singlet under the SM gauge group.The Higgs portal provides a unique gauge invariant and renormalizable interaction that couples S to the SM without the need of additional mediator particles.The model has first been considered in Ref. [1] and has since been extensively studied in the literature, see, e.g., Refs.[2][3][4][5][6][7][8][9].
The BSM part of the model Lagrangian reads (1) where H is the SM Higgs doublet and µ S , λ S , and λ HS are free parameters of the theory.In Eq. ( 1), the terms refer to (from left to right): the S kinetic term, the bare S mass, the S quartic self-coupling and the Higgs-portal coupling.The Lagrangian respects a discrete Z 2 symmetry that guarantees the stability of S.Under this symmetry all SM particles are assumed to be even while S is odd (S → −S).
After spontaneous electroweak symmetry breaking (EWSB), the Higgs boson acquires a vacuum expectation value, ⟨H⟩ = v 0 / √ 2, with v 0 = 246.2GeV, and Eq. ( 1) becomes: where h is the Higgs field and m 2 S = µ 2 s + 1/2 v 2 0 λ HS is the mass of S in the broken phase.
The quartic scalar self coupling term, proportional to λ S , is of importance for the stability of the electroweak vacuum and the perturbativity of the model [69] but does not affect the DM phenomenology.The latter is solely governed by the last two interaction terms in Eq. (2).
1. Tree-level annihilation Feynman diagrams of the SHP model.The symbol f refers to charged fermions, while the V refers to either W ± , Z.The central diagram includes also the u-channel.h q q q S S g g FIG. 2. Feynman diagram for the loop-induced annihilation SS → gg.The symbol q stands for quarks.
Thanks to these terms, S can annihilate into all SM particles through the Higgs portal with a coupling that is proportional to λ HS .The annihilation process is relevant for thermalization and freeze-out of S in the early Universe and can lead to the production of SM particles in our Galaxy today.Moreover, the Higgs boson may decay into the scalar S producing an observable signature at colliders.Finally, the scattering of S off quarks, mediated by h, could produce recoil events in direct detection experiments.
The SHP model is very simple, with only two parameters that are relevant for the DM physics: the physical DM mass m S and the Higgs portal coupling λ HS .This allows us to straight-forwardly constrain the parameters of the model through the relic density constraint and from direct and collider searches and astrophysical observations.In Figs. 1 and 2, we show the Feynman diagrams relevant for the annihilation process within the SHP model.The model contains one diagram for each fermion, one for each gauge boson W ± and Z and three diagrams for the annihilation into the Higgs boson.We do not show the channels with the production of γγ and γZ or γh that are loop induced, thus providing a subdominant contribution to the cross section.

B. Annihilation cross-section
Except for the channel with hh final states, the cross section σ multiplied by the Möller velocity v of DM pairs annihilating into SM particles i can be expresses as [70]: where Γ h→i ( √ s) is the partial decay width into state i of the SM Higgs boson evaluated at energy √ s.D h (s) is the Higgs propagator defined as: where Γ h is the total Higgs width, which includes all the kinematically allowed partial decay widths, as well as the invisible Higgs width Γ h,inv .For the former, we adopt the theoretical prediction from Ref. [71], Γ h,SM = 4.1 MeV.
For the latter, we employ the tree level result, given by: where m h is the Higgs mass.We use m h = 125 GeV as recently measured by ATLAS [72].The expression for σv written in Eq. ( 3) is particularly convenient because very precise theoretical calculations and measurements are present for the quantity Γ h→i (see, e.g., [71]).In Eq. ( 4), we write down the propagator within the commonly used fixed-width prescription.This is typically a good approximation except when Γ h is a rapidly varying function of s.This can, in particular, happen in the resonant region, m S ≈ m h /2, where the invisible decay channel opens up close to the resonance [73].In this case, the running of the Higgs width has to be taken into account by replacing m h in Eq. ( 5) with the centerof-mass energy √ s of the process.However, this effect is only relevant for sizeable coupling, as quantified at the end of Sec.III.
For the indirect detection signals, the cross section is averaged over the DM particle velocity distribution, for which we adopt the Standard Halo Model (see, e.g., [74]).However, the cross sections in the SHP are essentially s-wave, i.e. they do not depend on the velocity.The only region of the parameter space where the velocityaveraged annihilation cross-section ⟨σv⟩ depends on v is very close to the Higgs resonance.While this region is 10 5 10 4 10 3 10 2 10 1 10 0 10  90 GeV is dominated by the W ± channel while ⟨σv⟩ for mS = 330 GeV takes most of its contribution from the Higgs production.Right panel: Ratio between the annihilation cross-section for the different channels (⟨σv⟩i) and the total one (⟨σv⟩tot) as a function of the DM mass, computed for a fixed λHS = 0.01.
of high importance to our analysis, we nevertheless find that the velocity-dependence is relevant only for velocities significantly larger than the typical DM velocity in the Galactic halo.Therefore, for the typical Galactic velocities, which are of the order of 10 −3 c, the velocitydependent term at the resonance contributes only up to 1%, while being negligible otherwise.
The expression of the annihilation cross-section for all possible annihilation channels is provided in, e.g., Ref. [75].Here, we simply summarize the formulae for the relevant leading-order cross-sections and discuss their dependence on the parameters m S and λ HS .The scaling of the velocity averaged annihilation cross-section for different DM masses is shown in Fig. 3 and it will be used to guide the discussion.The annihilation cross-section into a pair of fermions is given by (see, e.g., [75]): where N c = 3 for quarks and N c = 1 for leptons, and m f is the fermion mass.If the S mass is above the Higgs resonance (m S > m h /2), Γ h is independent from λ HS and the annihilation cross-section into fermions scales as λ 2 HS (see Fig. 3 for m S = 90 GeV).This holds also for λ HS < 0.01 when m S < m h /2, because the decay of the Higgs boson into DM particles is negligible with respect to the decay into SM particles, and Γ h remains approximately equal to the SM one.Instead, for m S < m h /2 and λ HS > 0.01, Γ h gets a contribution which increases proportional to λ 2 HS [see Eq. ( 5)], thus adding an additional dependence on λ HS in the cross section, so that σ ∝ λ 2 HS /(1 + kλ 4 HS ).Note, however, that the second term in parentheses is subleading unless λ HS approaches the non-perturbative regime in which the validity of the expression becomes questionable.
The cross section into W ± and Z gauge bosons is given by: where V stands for W ± and Z.Most of the contribution from the annihilation into W ± and Z comes from DM masses above the resonance for which Γ h does not depend on λ HS .Therefore, the annihilation cross-section into W ± and Z scales as λ 2 HS in the most relevant parameter space (see Fig. 3 for m S = 90 GeV that is dominated by the W ± channel).
The cross section for annihilation into a pair of Higgs bosons gives rise to a more lengthy expression that we omit.However, in the non-relativistic limit, s → 4m 2 S , it reduces to: .
(8) In this case, for λ HS < 0.1 the annihilation cross-section scales as λ 2 HS .In contrast, for larger couplings, the polynomial term in the numerator becomes more relevant than all of the other terms, and the cross section scales as λ 4  HS .This is shown in Fig. 3 for the case of m S = 330 GeV, that is dominated by the annihilation into the Higgs boson.
We also take into account the annihilation into the gluon final states gg (see Fig. 2).This is a loop-induced process, which is typically suppressed with respect to the tree-level annihilation channels.In fact, its contribution is at most 7% for m S ≈ 60 GeV.We do not take into account the annihilation channels into γγ, γZ and γh because their contribution is at most at the few per mille level.
While we have explicitly shown here some expressions for the annihilation cross-section, useful to guide the discussion, we nevertheless calculate ⟨σv⟩ for all the annihilation channels using MadDM 1 [76][77][78][79], a plugin of MG5 aMC [80,81].We utilise the UFO implementation [82,83] for the SHP model.MadDM is linked to Pythia for hadronization and showering.In this work, we use MadDM version 3.2 [79] on top of MG5 aMC version 2.9.9 (LTS) and Pythia version 8.3 [84].
For the annihilation channels into massive gauge or Higgs bosons, we include their off-shell contributions.This is particularly relevant for W ± and Z final states, while it has very minor effect in the case of the h bosons.To this end, we calculate the cross section for the fourbody diagrams, where two bosons V V ′ are produced offshell from DM annihilation, subsequently decaying into four fermions f : The relative contributions of the different channels, ⟨σv⟩ i /⟨σv⟩ tot , is shown in Fig. 3 as a function of the DM mass, with λ HS = 0.01.For very small S masses, the annihilation happens predominantly into cc quarks and τ + τ − leptons.For m S = (5 − 50) GeV instead, the main contribution comes from b b annihilation.Finally, for large DM masses the main channels are the one with gauge and Higgs bosons.The effect of the off-shell contribution of the gauge bosons is particularly relevant for the annihilation involving W ± bosons, which becomes the main channel even below the on-shell production threshold m S < m W .
For DM masses above 5 GeV, the contribution of τ + τ − annihilation is larger than cc.This may seem to be counter-intuitive since, from Eq. ( 6), the cross sections for the two channels are σ cc ∝ 3m2 c and σ τ + τ − ∝ m 2 τ .Accordingly, one may expect However, taking into account the running quark masses in our calculation, we find m c ≤ 1 GeV [85] for m S ≥ 3 GeV, which leads to σ cc < σ τ + τ − .

C. Source energy spectra
The energy spectra dN p /dE of particles such as photons and antiprotons, produced from DM annihilation is an important ingredient for the indirect detection of DM (see Sec. VI).In order to calculate these spectra we use the code MadDM, which first evaluates the annihilation cross-sections.Subsequently, MadDM calls Pythia to generate the spectra for the different annihilation channels, and determines the final energy spectra by averaging the results obtained for each channel by ⟨σv⟩ i /⟨σv⟩ tot .
1 https://launchpad.net/maddm We include in the computation all channels (tree level and loop induced) that we consider for the calculation of the annihilation cross-sections (see Sec. II B).The spectra we generate take into account final-state-radiation (FSR) and electroweak emission of weak gauge bosons off fermions, which are included by Pythia.With respect to the computation of the tables provided in Ref. [86] there are two important additions.First, we include the contribution of off-shell W ± and Z, which could be very relevant between 60-90 GeV, as seen in Fig. 3. Second, we take into account the polarization of W ± and Z bosons.This effect is particularly relevant for the spectra of antiprotons.For all the technical details about FSR, electroweak corrections and inclusion of W ± and Z polarization see Appendix A.
We compute the DM spectra for masses between 2 GeV and 20 TeV generating 10 6 events per point.For DM masses above 110 GeV (where the hh channel introduces a dependence of the relative contributions of to annihilation on λ SH ) we vary λ SH between 10 −2 and 10. 2 In Fig. 4, we show the results for the DM spectra of γ rays and p for DM masses between 3 GeV and 10 TeV.The energy spectra, reported in terms of dN /d log 10 (x) , where x = K/m S with K being the kinetic energy, have a peak at an energy value that shifts towards lower values when increasing the DM mass.The p spectrum is decreasing for log 10 (x) < −2.Instead, the γ-ray spectrum has a flattening at around log 10 (x) < [−6, −3] depending on the mass.This change of shape is due to the contribution of the electroweak corrections such as FSR.

III. RELIC DENSITY
The DM average density has recently been measured by the Planck experiment to be Ω DM h 2 = 0.120 with an uncertainty at the level of 1% [87].In general, the theoretical calculation of the DM relic abundance requires to solve the Boltzmann equations for the DM phase-space density f S (p) in an expanding Friedmann-Robertson-Lemaître-Walker Universe [88,89]: where E and p are the energy and momentum of the S particle, and H is the Hubble expansion rate.The term C is the collision operator that takes into account all the interactions between DM and SM particles.In particular, C contains the operator for elastic scattering (C el ) and annihilation (C ann ).Elastic collisions are responsible for maintaining the kinetic equilibrium while inelastic collisions keep the chemical equilibrium.See Ref. [67] for a complete description of C el and C ann .For the canonical freeze-out of weakly interacting massive particles, a number of approximating assumptions are usually applied that significantly simplify Eq. ( 10) and its numerical solution.Most significantly, this concerns the assumption of kinetic equilibrium between DM and the SM bath throughout the entire process of chemical decoupling.The assumption that the distribution of S particles remains proportional to the thermal one, f S ∝ f S,eq , simplifies the partial differential equation Eq. ( 10) into an ordinary differential equation for the DM number density, i.e., the well-known Zeldovich-Okun-Pikelner-Lee-Weinberg equation [90,91]: where n S = g χ d 3 p/(2π) 3 f S (p) and ⟨σv⟩ T is the thermally averaged cross-section at temperature T .In the non-relativistic regime, the S particle phase-space density follows Maxwell-Boltzmann statistics f S,eq ≃ exp (−E(p)/T ) and ⟨σv⟩ T reads: where K i are the modified Bessel functions of order i.
The assumption of kinetic equilibrium is often well justified as elastic scattering processes can easily be orders of magnitude larger than the annihilation processes that initiate chemical equilibrium.This is due to the large number density of light SM particles DM can scatter off (compared to the Boltzmann suppressed number density of heavy DM particles during freeze out).However, this reasoning does not carry over to annihilation via a resonant Higgs in s-channel.First, elastic scattering processes are not resonantly enhanced and, hence, suppressed compared to annihilation.Second, as the Higgs-SM couplings are proportional to the SM particles masses, coupling to light particles -with unsuppressed number densities -are small.As a result, kinetic equilibrium can break down during chemical freeze-out and the above assumption is unjustified [67].In this case, the full version of the Boltzmann equation must be solved, with the inclusion of the elastic collision term.We compute the relic density for S, hereafter called Ω S h 2 , in both setups, labelled as: • fBE: the full solution of Eq. ( 10); • nBE: the solution of Eq. ( 11), i.e., assuming kinetic equilibrium during freeze-out.
In Ref. [67] the authors found that the relic abundance predicted following the approximated approach nBE can differ by up to an order of magnitude from the calculation obtained with the full Boltzmann equation, fBE.We performed the relic density calculation using the code DRAKE [68], developed by the authors of Ref. [67].This code solves the Boltzmann equation with both the models nBE and fBE.For a detailed explanation, see Ref. [68].
For the values of m S we considered, the freeze-out takes place at temperatures around a few GeV.That means we are close to the region in which the QCD phase transition takes place, after which quarks become confined in baryons and mesons.To take into account this effect, we considered two different physics models, which are implemented in the DRAKE code and studied in Ref. [67].These two models are extreme scenarios that are intended to bracket the uncertainties: The first model, named QCD A , represents the case maximizing the elastic scattering for which all quarks are free and present in the plasma (according to their equilibrium abundance) down to temperatures of T c = 154 MeV.The second model, named (QCD B ), minimizes the elastic scattering: only light quarks (u, d, s) can contribute, and  only considering temperatures above 4T c ≈ 600 MeV.This threshold ensures that hadronization effects are negligible.
For comparison, we also compute Ω S h 2 in the setup nBE with micrOMEGAs 3 [92][93][94][95].We consider the SHP model implementation shipped with the code.Both DRAKE and micrOMEGAs take into account the offshell contributions of W ± and Z gauge bosons, and the contribution coming from loop-induced annihilation into gluons.
The left panel of Fig. 5 shows the combination of parameters λ HS and m S for which the relic density Ω S h 2 corresponds to the measured value Ω DM h 2 = 0.120.We show the computation of the fBE model in the QCD B scenario with DRAKE.The first consideration we can do regards the resonance region: for m S ≈ m h /2 the value of λ HS decreases very quickly.This can be explained with the resonant enhancement of the annihilation crosssection for this particular parameter choice.As a result, the correct value of the relic abundance necessarily requires a very small value of the coupling λ HS .The right panel of Fig. 5 shows the parameter λ HS accounting for the correct relic density value as a function of m S in the resonance region.We show the computations performed with micrOMEGAs (which includes only the nBE model) and with DRAKE (for both the fBE QCD A and QCD B models).The difference between the fBE case and the normal approach nBE varies within a factor 0.5 and 2, and it is present in the range of DM masses between (50 − 67) GeV.In particular, around 57 GeV, the 3 https://lapth.cnrs.fr/micromegas/fBE result is approximately a factor of 2 larger than the nBE one, while on the Higgs pole it differs by a factor of ∼ 0.5.In the figure we also show the difference between the QCD A and QCD B assumptions, which is maximal in the range (55 − 60) GeV, where it is at the level of 40 − 50%.This region represents the uncertainty due to the QCD phase transition.This is our main result for the computation of the relic density, which properly takes into account the kinetic and thermal parts of the Boltzmann equation, thanks to the implementation of the fBE model, additionally accounting for the QCD phase transition effect in the resonance region.
There are further corrections to the annihilation channel that can potentially become important close to the resonance.As we will discuss in Sec.IV in the context of LHC constraints, for m S very close to m h /2, the centreof-mass energies relevant in the annihilation process are around the SS threshold.Accordingly, the change of the Higgs width due to the opening of the SS decay channel can become relevant and the fixed-width calculation can break down calling out for a running-width description [73].However, for annihilation, we find that this effect is larger than a few percent only for λ HS > 0.3 in the region |m S − m h /2| ≈ 0.1 GeV.Accordingly, in contrast to its relevance for LHC constraints, the effect can be neglected in the computation of the relic density since considerably smaller couplings are required in the resonant regime.This is true even for the scenarios in which S constitutes a subdominant fraction of DM considered in Sec.VII.We neglect thermal corrections to the Higgs width and virtual corrections to the annihilation into an on-shell Higgs which have recently been computed in Ref. [96].However, these computations assume kinetic equilibrium.
Finally, we mention an alternative mechanism to generate the DM density within the model.For very small (or feeble) Higgs-portal couplings, DM may never thermalize with the SM bath.In this case, DM can be produced via freeze-in [97][98][99].Explaining the measured relic density in this scenario requires λ HS ∼ 10 −12 −10 −11 for m S = (50 − 70) GeV.Due to the non-efficient annihilation processes in this regime, an initial abundance from processes taking place during reheating cannot be diluted as it is the case for thermalized DM discussed above.Hence, this scenario faces further constraints from physics of the very early Universe.Furthermore, the annihilation rate is way below the one providing a detectable signal in astroparticle data, such as the GCE (see Sec. VII).We do not consider this case here.

IV. COLLIDER SEARCHES
Being independent of astrophysical uncertainties and cosmological assumptions, collider searches for invisible decays of the SM Higgs boson into SS constitute a robust probe of the model.For m S < m h /2, the singlet scalar S can contribute to the invisible width Γ h,inv of the Higgs boson and we can directly reinterpret the limits derived by the experimental collaborations.The invisible branching ratio is: The latest 95% CL upper limits on B h,inv are 0.13 [100] (ATLAS, preliminary) and 0.19 [101] (CMS, published).
In this work, we are interested in the resonant region, m S ≈ m h /2, where contributions of an off-shell Higgs boson production are relevant.Furthermore, as pointed out in Ref. [73], the total Higgs boson width can become a rapidly varying function of the invariant mass due to the opening of the Higgs boson decay channel into SS requiring a computation beyond the fixed-width description.This calls out for a more careful reinterpretation of the experimental results.To this end, we follow the method presented in Ref. [73].
Assuming the factorization of the Higgs boson production and decay channels, in an analogue way to Eq. ( 3), we can express the cross section for DM pair production at the LHC as where σ h (s) is the Higgs boson production cross-section for an (off-shell) Higgs boson with invariant mass √ s.This equals the on-shell production cross-section of a SM-Higgs-like scalar ϕ with mass m ϕ , σ h (s = m 2 ϕ ) = σ ϕ (m ϕ ).Ref. [100] reports a 95% CL upper limit on the product of the cross section times the invisible branching ratio for such a scalar.Using this result, we can employ Eq. ( 14) 6. 95% CL upper limits for the parameter λHS as a function of the DM mass mS.We show the constraints obtained from a reinterpretation of the ATLAS search in Ref. [100] and the projected limits for the HL-LHC and HE-LHC at 14 and 27 TeV, respectively, taken from [73].
to obtain limits on the DM coupling in our model without the use of Monte Carlo simulations.We use the running width in the Higgs propagator, i.e. replacing m h → √ s in Eq. ( 5) (see the appendix of Ref. [73] for further details).We take the theoretical predictions for the cross section σ ϕ (m ϕ ) from Refs.[71,102].
The result for the 95% confidence level (CL) upper limit on λ HS is shown in Fig. 6 with a black, dashed curve, as a function of m S .For m S well below m h /2 this result resembles the quoted upper limit on B h,inv of 0.13, obtained in [100], while it deviates strongly from the so-obtained limit for m S ≳ m h /2.Notice that the constraint is about a factor of 2.5 stronger than the one derived in Ref. [73] on the basis of Ref. [101].We also show the potential for future improvement of the limit at the HL-LHC with 14 TeV and a HE-LHC upgrade with 27 TeV center-of-mass energy taken from Ref. [73].

V. DIRECT DETECTION SEARCHES
Direct detection experiments such as LZ and XENONnT can provide very tight constraints on the cross sections for the interactions of DM particle with nucleons [103].Direct detection experiments typically measure upper limits on the elastic scattering cross-section off nucleons, as a function of the DM mass.For the SHP model, the cross section is spin-independent only.Both the XENONnT and the LZ experiments have recently released the tightest constraints so far on the spinindependent cross-section [62,63] given as a 90% CL upper limit on σ SI which are at the level of 10 −47 cm 2 for DM masses of the order of few tens of GeV.The upper limits have been derived assuming the Standard Halo Model.The local DM density ρ ⊙ is conventionally fixed at 0.3 GeV/cm 3 , while DM velocity follows a Maxwellian velocity distribution, with velocity dispersion 220/ √ 2 km/s, escape velocity v esc = 544 km/s, and Earth velocity of v E = 232 km/s.XENONnT and LZ are expected to reach upper limits at the level of 10 −48 cm 2 [103].
The spin-independent cross-section in the framework of the SHP model is (see [70]): where m N is the nucleon mass.
The term f N parametrizes the Higgs-nucleon interaction: where the sum is made over all quark flavours q with mass m q .A typical value of f N ≈ 0.3 is often adopted [70].
To compute the spin-independent cross-section for the SHP model, we used two numerical codes: MadDM [77] and micrOMEGAs.The results found with mi-crOMEGAs are larger than the one of MadDM by a factor of 5%, independent from the DM mass.The main features responsible for the difference are associated to the QCD corrections and to the contribution of higher twist operators, which are not taken into account in MadDM 4 .
In Fig. 7 we show the upper limits found using the available data from the LZ [62] and XENONnT [63] experiments.In particular, LZ can probe values of the coupling parameter λ HS as low as about 10 −3 for DM masses of a few tens of GeV.However, the local DM density profile and its velocity distribution is characterized by sizable uncertainties.For DM masses below 10 GeV, the uncertainties on the velocity distribution play an important role, while the ones on the value ρ ⊙ are the dominant source of uncertainty at higher DM masses (see [104]).Since the parameter space we are considering consists of DM masses above 10 GeV, we take into account only the uncertainties on ρ ⊙ , and we consider it to vary in the range (0.2 -0.6) GeV/cm 3 [105][106][107][108].In Fig. 7, left panel, we show the uncertainty band associated to ρ ⊙ .In the same plot we also show the projected sensitivity for the future experiment DARWIN [65], which will reduce the upper limits by nearly a factor of 4 for the considered mass range.
Assuming that the total DM relic density can be associated to only one single DM candidate, in our case the S scalar, the differential rate of detection dR/dE is proportional to σ SI SN ρ ⊙ /m S , where ρ ⊙ is the local DM mass density.On the contrary, assuming that the candidate S constitutes only a fraction of the full DM component, i.e. that the relic density of the particle S is a fraction of the total DM relic density, we have to apply an appropriate rescaling to the cross section upper limit.Indeed, in this case, the local DM density of the DM candidate S is lower, and can be simply assumed to scale as: Assuming that there is no difference in how the particle S and the other DM components cluster in the Galaxy, the local energy density of S can be written as As a consequence, the upper limits on σ SI associated to S can be computed as the parameter space of m S , λ HS , satisfying the condition: where σ UL EXP is the 90% CL upper limit of the experiment under consideration.In Fig. 7, right panel, we show the results obtained by including the effects of considering a fraction of DM into the calculation of the relic density.Using the LZ upper limit, the only region allowed corresponds to S masses in the range (57−62.51)GeV.Additionally, we also consider the projected DAR-WIN upper limits and the remaining narrow region is m S = (61 GeV − m h /2), i.e below m h /2.Over the resonance region, the couplings compatible with the direct detection constraints are much larger than what is shown in Fig. 7.For example, for m S = 60 GeV, the upper limits on the coupling parameter increases by a factor of about 100.Therefore, the requirement of achieving the correct relic density for the S particles in the calculation of the direct detection has the consequence of restricting the mass range compatible with the data upper limits but it provides weaker constraints on the λ HS parameter in the same mass range.Future experiments seem to not be able to rule out values of λ HS smaller than 0.1 if the DM mass remains slightly below m h /2.

VI. INDIRECT DETECTION
Current experiments measuring cosmic particles, such as γ rays, neutrinos or CRs, count the indirect detection of a possible DM signal among their most important science cases (see, e.g., [42] for Fermi-LAT).In particular, these experiments focus their searches on the detection of fluxes of the rarest astrophysical particles such as: photons (from radio to γ rays), neutrinos and antineutrinos ν, antiprotons p, positrons e + and antinuclei.Among these particles, photons and, to an even S h 2 = 0.12 FIG. 7. Left panel: The 90% CL upper limits for the parameter λHS as a function of the DM mass mS obtained with direct detection experiments.We recast the exclusion limits of LZ [62] and XENON1T [63] along with the projected sensitivity for the future experiments DARWIN [65] for the SHP model.The grey region represents the uncertainty in the upper limit induced by the uncertainty on the local DM density ρ⊙, which we vary in the range (0.2 -0.6) GeV/cm 3 .Right panel: Same as left but with the constraints re-scaled on the basis of the actual value of ΩSh 2 following the assumption that the DM might be subdominant or overabundant, i.e., by using Eq. ( 19).We show the allowed region obtained using the upper limits from the LZ experiment [62] (green region) and the projected sensitivity for the future DARWIN experiment [65] (orange region).We also show the combination of λHS and mS that provides the correct relic abundance and we display the region where the particle S is overabundant with a grey region.
larger extent, ν have the great advantage to travel almost undisturbed across the Universe and are able to provide directional information.This allows to focus the DM search in the direction of astrophysical sources which are expected to have very large density of DM.Among the most promising targets, we mention the Milky Way's Galactic center, dSphs and the clusters of galaxies.Current experiments can detect astrophysical neutrinos for energies above 1 TeV and with quite low statistics (see, e.g., [109]).These two aspects make them not suitable for our scopes.Concerning photons, we focus on the γ rays detected by Fermi-LAT and we use the data available for the GCE from Refs.[27,28] and the recent analysis of DM searches in dSphs from Refs.[29,51].Finally, we consider among the CRs only antiprotons because the SHP predicts a large hadronic production of these particles, while the e ± production is less important.Finally, we won't consider cosmic antinuclei because currently no firm detection has been published so far and upper limits are much weaker than other probes (see, e.g., [110]).
A. The γ rays from the Galactic center and dwarf galaxies

Model
The γ-ray emission from DM particle interactions includes two components: the direct and indirect production.The first one is also called prompt emission and it is due to the direct production of γ rays through an intermediate annihilation channel.The prompt emission is theoretically calculated as follows: where ρ ⊙ is the local DM density, and r ⊙ is the distance of the Earth from the center of the Galaxy.We use r ⊙ = 8.12 kpc, as measured recently in Ref. [111].The term J is the geometrical factor averaged over the solid angle ∆Ω spanned by the region of interest (ROI) considered in the analysis.This quantity is calculated as the integral performed along the line of sight (l.o.s.) l of the squared DM density distribution ρ 2 normalized by ∆Ω: The term (dN γ /dE) f is the γ-ray spectrum from DM annihilation for a specific annihilation channel labeled as f and ⟨σv⟩ f /⟨σv⟩ is its relative contribution of the annihilation cross section into a specific channel f with respect to the total annihilation cross section.As explained in Sec.II C we calculate the source spectra in the SHP model using MadDM.
In case of DM particles annihilating into leptonic channels, i.e. e + e − , µ + µ − and τ + τ − , there is a secondary production of γ rays that could be relevant.This involves e ± produced from the prompt emission that can subsequently generate γ rays through inverse Compton scat- tering on the interstellar radiation fields photons.This component is particularly relevant for the Galactic center where the density of the starlight and dust components of the interstellar radiation fields are roughly a factor of 10 higher than their local density (see, e.g., [112]).As demonstrated in Ref. [29], hadronic annihilation channels have a contribution from inverse Compton scattering relevant only for energies below 1 GeV, for m S between tens of GeV up to 100 GeV, which is the most interesting region in our analysis.For these energies, the uncertainties in the GCE flux are very large and the significance of detection is small.Therefore, we decide to neglect the contribution of secondary γ rays coming from DM annihilation.
The exact DM density distribution in the center of the Galaxy is not well known.In Ref. [29], the authors considered an hybrid approach to estimate the uncertainty on the density profile.They took into account two factors: the data on Galaxy rotation curve, which constrains ρ beyond a few kpc from the Galactic center; and the GCE data, which fixes the density shape between about 0.1 and few kpc under the assumption that the GCE is originated by DM annihilation.Using this strategy they identified three models, labeled as MIN, MED MAX, that bracket the uncertainty on the geometrical factor.In particular, they have found that the variation for the J value between the MIN and MAX models is about a factor of 7. We considered this uncertainty in our results.Instead, the uncertainty on J for dSphs is included directly in the analysis of the γ-ray data from these objects.In fact, this parameter is treated as a nuisance parameter of the likelihood fit, using a gaussian prior with an average taken as the observed geometrical factor and with the 1σ error as the width (see [29,51] for more details).

Results
We first provide the results of the fit to the GCE data found in [27,28].We perform a fit to the data taking into account statistical and systematic errors.The latter include the uncertainties due to the choice of the Galactic interstellar emission model.The DM mass m S and the annihilation cross section ⟨σv⟩ are free parameters of the fit.The left panel of Fig. 8 shows the result of the fit to the data found in [28], which gives the best-fit values for the DM parameters: m S = (58 ± 6) GeV and ⟨σv⟩ = (2.4 ± 0.4) × 10 −26 cm 3 /s.The best-fit value for the coupling parameter is between 0.03-0.18.It varies significantly because the best-fit masses are close to the Higgs resonance for which the ⟨σv⟩ increases dramatically, and the required coupling decreases.This is visible in the right panel of Fig. 8, where we show the χ 2 profile as a function of the parameters m S and λ HS .When m S ≈ 50 GeV, the value of λ HS that fits the data is approximately (1−2) × 10 −2 .Close to the resonance m h /2, the best-fit of the coupling decreases to values of the order (1−2) × 10 −4 .This is due to the fact that by fixing λ HS , the ⟨σv⟩ increases significantly when m S ≈ m h /2.For masses around (70−80) GeV, the best-fit value of the coupling increases to values of the order 10 −1 .The χ 2 profile we find using the data released in the other reference considered in this paper for the GCE, i.e.Ref. [27], is similar to the one shown in Fig. 8 and obtained from the GCE SED obtained in Ref. [28].As a final remark, we remind that the geometrical factor for the Galactic center is not well known.As seen before, in Ref. [29] an uncertainty of about a factor of 7 in the J value has been estimated between the models describing the GCE density profile.This translates into a systematic also in the best-fit of λ HS .In particular, the calculation of the γ-ray flux is proportional to ⟨σv⟩ • J .That means a variation of J can be reabsorbed by an opposite variation of the annihilation cross section, namely in a variation of the value of the coupling λ HS , where we have ⟨σv⟩ ∝ λ 2 HS .Therefore, a systematic of a factor of 7 in the geometrical factor causes a systematic in the value of λ HS of about √ 7 ≈ 2.6.We also analyzed the γ-ray constraints coming from dSphs.In particular we use the results of Refs.[29,51].In the former, a sample of 48 dwarf galaxies, both classical and ultrafaint, have been analyzed using 11 years of data and assuming that the objects are pointlike (i.e.smaller than the Fermi-LAT point-spread function).In the latter, 12 years of data in the direction of 22 dSphs have been analyzed considering a spatial extended templates for the DM flux.In both cases no clear detection has been found with a maximum value of the test statistic of the order of 9. Using the data released in Refs.[29,51] we determined the χ 2 profile as a function of the m S and ⟨σv⟩ as we did for the GCE.The peak of the signal is located at masses around (50−300) GeV and ⟨σv⟩ ≈ (0.4−6) × 10 −26 cm 3 /s.

B. Antiproton flux
DM annihilation can induce a primary flux of p.The source term which denotes the differential production rate of p per volume, time and energy reads as: where ( dN p/dE ) i is the spectrum of p produced from DM particle interactions and it depends on the specific annihilation channel assumed and labeled with i.In addition, there is an astrophysical antimatter background which originates from the scattering of CR protons and nuclei on the interstellar matter.
To properly calculate the flux of secondary and DM antiprotons one should calibrate the galactic propagation parameters such that the spectra of primary particles reproduce AMS-02 observations.We use the results presented in Ref. [113] where the authors have first fitted the data sets of proton, helium, and the antiproton-toproton ratio from AMS-02 collected over 7 years [114] and complement these data sets with low-energy data for protons and helium from the Voyager satellite [115].Their propagation model includes diffusion (parametrized by a smoothly broken power-law in terms of the particle rigidity), reacceleration, energy losses, secondary production/fragmentation processes, and solar modulation modelled is a force field approximation (see Ref. [113] for details).By adding a possible DM contribution to the antiproton source spectra, our analysis provides a χ 2 profile as a function of the parameters m S and λ HS of the SHP model.We combine these results for p, published in Ref. [113], with the ones found for γ rays from the Galactic center and dSphs.The authors recently followed up with an updated analysis [116].In that work, they have tested two models, labelled as INJ.BRK and DIFF.BRK.The former describes the injection spectra of the primary CRs with a broken power law using different slopes for p and He.The diffusion coefficient is modeled as a broken power law with a single break.The latter model, instead, utilizes a single power law for the spectra, while the diffusion coefficient has both a low and an high-energy break.In this analysis, the authors perform fits to the AMS-02 antiproton data with and without including the possible correlations [59] in the data.The INJ.BRK model provides evidence for a possible DM contribution similarly to what has been found in Refs.[110,113,117].
For completeness, we considered both models in our analysis.In particular, the INJ.BRK model results in very minor changes on our constraints, while the DIFF.BRK model shows that there is no evidence for a DM contribution for masses below 100 GeV.As a consequence, the latter model allows us to impose tight constraints on the SHP model parameter space.We will discuss this possibility in Sec.VII C.

VII. COMBINED RESULTS FROM DIRECT, INDIRECT, COLLIDER SEARCHES AND COSMOLOGY
A. The S particle as dominant dark matter component In this section, we combine the results found previously using γ-ray and p cosmic data, collider, direct detection constraints and cosmological measurements.We first combine the results obtained using γ-ray and p flux data by summing the χ 2 profile defined as a function of the parameters m S and λ HS .In Fig. 9 we show the best-fit region from astroparticle data compared with the constraints from collider, direct detection and relic density.For the direct detection we show the upper limits on λ HS obtained with the LZ data (see Sec. V) while for collider searches we use the upper limits found with CMS data for the vector boson fusion (see Sec. IV).In this figure as well as in Fig. 10, we assume that the S particles make up 100% of the DM relic density.This assumption enters the calculation of the indirect and direct detection signals whereas the green band shows the coupling that provides the respective relic density through thermal freeze-out assuming a standard cosmological history.There are two intersections of the green and color-shaded bands indicating compatibility with the relic density and astroparticle data, respectively.The first one is at a DM mass m S ≈ m h /2 and λ HS ≈ (1.4 − 1.7) × 10 −4 .This is compatible with the collider and direct detection constraints.The second one is for m S = (63 − 69) GeV and λ HS = (2 − 6) × 10 −2 .However, this part of the parameter space is ruled out by direct detection constraints.Therefore, the only region of the parameter space that fits well cosmic data, compatible with relic density and  consistent with collider and direct detection constraints is for m S ≈ m h /2, to be precise (10−20) MeV smaller than m h /2, and λ = (1.4−1.7)× 10 −4 .In this region of the parameter space the GCE is fitted with a reduced χ 2 of about 0.8, which implies that the fit is statistically appropriate.
As explained in Sec.VI A, there is a quite large uncertainty related to the exact DM density in the center of the Milky Way, which in Ref. [29] has been estimated to be included among the MIN, MED and MAX models.In Fig. 10, we show the combined results obtained with the MIN and MAX models observing that a change of the DM density in the center of the Galaxy has little to no effect on our conclusions.With these models, there is still a region very close to the Higgs resonance which fits well cosmic fluxes data and which is compatible with the collider and direct detection experiments.The only quantity that changes is the value of λ HS that fits the cosmic fluxes data and for which we have the correct relic abundance.By considering the MIN and MAX model the coupling can vary in the range (1.2−2.0)× 10   11.The same as in Fig. 9 but for various values of the ξ parameter, as labelled in the different frames.

B. The S particle as subdominant dark matter component
In this section, we discuss the possibility that the S particle is not the dominant component of DM density.We therefore allow paramenter-space configurations for which Ω S h 2 accounts from 1% up 100% of the measured dark matter average density Ω DM h 2 (i.e.0.01 ≤ ξ ≤ 1).We verify whether the relic density regions are still compatible with astroparticle data best-fit regions and with direct detection and collider constraints.In Fig. 5, we show the parameter space of the SHP model providing a fraction of the DM relic density.In order to find the best-fit region for λ HS and m S that fits cosmic fluxes data, we consider the following.If the S particle constitutes a fraction ξ of the DM relic abundance, according to Eq. ( 17), the DM Galactic density should be rescaled by ξ.The geometrical factor J is proportional to ρ 2 .This implies that J ∝ ξ 2 .Since Φ ∝ J ⟨σv⟩ ∝ ξ 2 ⟨σv⟩, to obtain the same flux, ⟨σv⟩ has to increase for decreasing ξ according to ⟨σv⟩ ∝ ξ −2 .Since ⟨σv⟩, in turn, scales with λ 2 HS , for λ HS < 0.1 (see Fig. 3), we can just rescale the values for λ HS for the combined fit to the cosmic particle flux data shown in the previous subsection by 1/ξ.
In Fig. 11, we show the results for the combined analysis for ξ = 0.3, 0.1, 0.03, and 0.01.When ξ = 0.3 (S makes up 30% of total DM density) we see that the bestfit region to cosmic data matching the relic abundance corresponds to m S ≈ m h /2 and λ HS = (2.3−2.7)× 10 −4 , which is consistent with the direct detection upper limits.The second region with mass around m h /2 + 0.20 GeV and λ HS = (0.8 − 1.8) × 10 −2 is instead excluded by LZ constraints.For ξ smaller than about 20%, also the higher mass region of compatibility between cosmic data and relic abundance becomes consistent with direct detection upper limits.This is visible in Fig. 11 for the case with ξ = 0.1 where there is still a region of compatibility at about m S ≈ m h /2 and λ HS ≈ 3 × 10 −4 and the second one which is at m S ≈ m h /2 + 0.05 GeV has λ HS ≈ (0.5−1.2) × 10 −2 .The smaller the value of ξ, the closer the high mass compatibility region is to the right side of m S ≈ m h /2 with at the same time smaller values of λ HS .For ξ < 0.02 the two regions merge.This can be  FIG.13.Best-fit regions of the combined cosmic particle flux data compatible with thermal freeze-out (corresponding to the intersections of the green and color-shaded bands in the plots of Fig. 11) in the ξ-λHS plane.The green (blue) bands denote the regions where MS is slightly smaller (larger) than m h /2.We also report the upper limits from direct detection covered by LZ and the future prospects for DARWIN.
seen in Fig. 13 showing the two best-fit regions of the cosmic particle flux data for λ HS as a function of ξ demanding compatibility with thermal freeze-out (corresponding to the intersections of the green and color-shaded bands in the plots of Fig. 11).For values of ξ smaller than 1% there is no region of the parameter space for which the best-fit from cosmic data is compatible with the relic density calculation.Therefore, if ξ ≥ 0.20 the only region fitting the cosmic data and compatible with the relic density and direct detection constraints is for m S ≈ m h /2 with values of λ HS of the order of (1.5−4) × 10 −4 .Instead, if ξ = 0.01 − 0.20, also slightly larger m S values are allowed, of the order of 0.05 GeV larger, while the coupling can be as large as 10 −2 .
To summarise our results, in Fig. 12, we have plotted the best-fit region and constraints while rescaling all astrophysical signals according to the value of ξ obtained from the relic density computation for each point in the parameter space.In the green shaded region, Ω S h 2 > 0.12, i.e. this region would over-close the Universe and it is hence excluded.Above the shaded region, ξ < 1 with ξ decreasing towards larger λ HS .The orange shaded region is excluded by direct detection constraint from LZ.As in Fig. 9, the grey dashed line denotes the (ξ-independent) constraints from the LHC with the region above the line being excluded.In the right panel of Fig. 12, we zoom into the resonant region revealing the narrow parameter space preferred by the observations.The projected sensitivity of DARWIN is shown in Figs.11 and 13 as the purple dot-dashed curves.While the region with m S ≲ m h /2 and smaller coupling remains out-of-reach for direct detection experiments in the foreseeable future, the region with m S ≳ m h /2 (blue shaded area) can be probed with the sensitivity of DARWIN.
The results presented in this section are consistent with the ones reported in Ref. [24].

C. Results without the GCE and upper limits from cosmic antiprotons
An alternative interpretation of the GCE states that it originates from the cumulative flux coming from a population of millisecond pulsars located in the Galactic bulge [35][36][37][38].If this is the case, tight constraints on a possible DM contribution in the Galactic center can be found (see, e.g., [118]).Additionally, the analysis of a possible DM contribution to the cosmic ray flux of antiprotons have not found any signal, and upper limits on the annihilation cross section has been imposed in Ref. [116] with the model DIFF.BREAK.
In this section, we report the results of our analysis assuming that the DM is not responsible for the GCE.We take into account the results obtained by Ref. [116] with the model DIFF.BREAK.We can not use the results of Ref. [118] since they do not consider the SHP model.
In Fig. 14, we show the results obtained with a combined analysis of dSphs, as presented in Sec.VI A, and with antiprotons, following the results published in Ref. [116].We considered both DIFF.BREAK and INJ.BREAK models, with and without taking into account correlations between data points.In particular, we display the 95% CL upper limits for λ HS as a function of the DM mass obtained for the two propagation setups used in the antiproton analysis and the two cases of data correlations.The upper limits rule out the region of the parameter space with m S ≥ m h /2 because for those masses the values of λ HS which satisfy the relic density are above the upper limit curves.Instead, for m S < m h /2 the allowed region contains masses below m h /2−0.05GeV and λ HS smaller than a few times 10 −4 .

VIII. CONCLUSIONS
In this paper, we performed a detailed and updated phenomenological study of the SHP model combining data from indirect and direct detection, the LHC, and the DM density measurement.We pinpointed the parameter space that is both allowed by the constraints imposed by this data as well as preferred as an explanation of the long-standing Fermi-LAT GCE and the AMS-02 antiproton excess.The former set of data alone leaves the resonant region m S ∼ m h /2 as one of only two allowed windows of parameter space -the other one being the high-mass region, m S ≳ 3 TeV, which, however, requires large couplings λ HS ≳ 1 and will entirely be probed by future direct detection experiments, e.g., with DARWIN.In the resonant region, in contrast, very small couplings are sufficient to match the measured relic density containing a spot with maximal resonant enhancement that is out-of-reach for any direct detection and collider experiment in the foreseeable future.This leaves indirect detection as the only handle to conclusively probe this scenario.
Accordingly, we further narrow down the parameter space in the resonant region by testing its compatibility with the Fermi-LAT GCE and the AMS-02 antiproton excess independently pointing to a DM mass around m h /2.We utilize recent results from Refs.[27,28] and [113] based on 11 years of Fermi-LAT data and 7 years of AMS-02 data, respectively, thereby improving over previous results obtained within the model.Furthermore, special care has been taken in the precise computation of the relic density by considering effects of kinetic decoupling that are highly relevant close to the resonance.Two viable subregions fit both excesses while predicting the right relic density: one slightly below and one just above m h /2 among which the latter is, however, excluded by direct detection due to the larger coupling required.When allowing S to only constitute a fraction of the DM density -while still being considered exclusively responsible for the observed indirect detection excesses -more possibilities open up.In particular, due to the different scalings of indirect and direct detection rates with the DM fraction, the latter subregion becomes allowed for a fraction of 10% or less.Both subregions provide good fits down to a DM fraction of 1% below which the relic density and the indirect detection signals become incompatible.
Providing couplings in between 10 −2 and 10 −4 , both regions are extremely hard to probe with future experiments.We showed that a HL-LHC and even a 27 TeV HE-LHC could not test the model in the laboratory.However, with the sensitivity of DARWIN, the subregion with m S ≳ m h /2 and coupling λ HS ≳ 2 × 10 −3 can be probed entirely.
A Python package enabling fast calculations of the dark matter spectra, relic density, as well as direct detec-tion and collider constraints within the model is available on Github.
Appendix A: Dark Matter energy spectra calculated with MadDM In this section we provide the technical details for the calculation of the DM energy spectra with MadDM.In particular, we refer to the computation of the quantities dN p /dE with p being γ rays, positrons, antiprotons or neutrinos, that enter in the calculation of the flux of particles produced from DM annihilation (see Sec. VI).
As a first step, given a value of the DM mass m S and coupling λ HS , MadDM calculates the cross section for all the annihilation channels.Then, it obtains the relative contribution of each i-th channel as ⟨σv⟩ i /⟨σv⟩, and it calculates the final spectrum of a certain particle p as: summing over the different annihilation channels labeled by the index i.During the calculation of the annihilation cross sections, MadDM calls the MadEvent generator to simulate events of the hard process, including the kinematics of the produced particles, eventually saving all the information in a Les Houches Event (LHE) file [119].
Event numbers are generated proportional to the relative contributions of the channel i, ⟨σv⟩ i /⟨σv⟩ tot .
In the next step, MadDM internally calls the Pythia Monte Carlo generator, which takes as input the LHE file.In particular, Pythia reads the LHE files and performs a Monte-Carlo simulation of the hadronization and showering of the DM annihilation final states, starting from the kinematics defined in the LHE file itself, and producing final energy spectra of stable particles, such as photons, neutrinos, e ± and p.At this stage, loop-induced processes are handled differently: annihilation cross section calculation and related event generation, as well as the spectra simulation through Pythia, are performed independently and separately from tree-level processes.Eventually, the spectra obtained from both the tree-level and loop-induced processes are combined according to their respective branching ratios.
The energy spectra we generate take into account FSR, which is included by Pythia using the command PartonLevel:FSR=on, while we turn off initial-stateradiation (PartonLevel:ISR=off), which is not present for DM annihilating particles, and multi-particle interactions (PartonLevel:MPI=off).We turn on the emission of weak gauge bosons off fermions which is part of FSR by using TimeShower:weakShower = on.Our procedure differs from the one used in the Ref. [86] where the authors have generated DM spectra in a model independent framework.In particular, that approach implies the process is created through a fake resonance with energy equal to twice the DM mass, and making it decay into the channel of interest.Some important effects are not taken into account if following this method.The first is related to off-shell contributions from the massive gauge bosons.As we have seen in Fig. 3, the contribution of off-shell W ± and Z could be very relevant between  GeV.The second effect is related to the polarization of W ± and Z bosons.By using the resonant mechanism as in Ref. [86] the spin 1 gauge bosons are unpolarized, resulting in the final fermions produced after their decay to acquire the inaccurate kinematics.Instead, in our calculation, we include the polarization of the gauge bosons because, thanks to MadDM, we generate events up to the fermions produced from the decay of the W ± and Z particles.When Pythia reads the LHE files, the kinematics of the final fermions inherited the fact that the W ± and Z have a polarization.
We show in Fig. 15 the effect of the polarization of the W ± and Z bosons for the production of antiprotons and γ rays.In particular we consider three DM masses: 100, 300 GeV and 3 TeV.The effect is more relevant for higher DM masses in the p channel, with respect to γ rays.In fact, for photons and m S = 100 GeV the effect accounts for a few % correction up to energies very close to the DM mass where the unpolarized spectra are larger by a factor between (20−25)%.Increasing m S , the effect on the γ-ray spectra starts to be visible also at lower energies where the differences can be of the order of (10−20)%.Instead, for the production of p, the effect could be larger and reaches even (40−50)% at small energies and large values of m S .In any case, the effect is small at the peak of the spectra, in units of x 2 ( dN /dx ), where x = E/m S , and increases in the low and high-energy tails.

FIG. 4 .
FIG. 4. Left panel: Energy spectra ( dN /d log 10 (x) ) of p for the SHP model plotted as a function of x = K/mS, where K is the kinetic energy.We show the energy spectra for different DM mass values ranging from 3 GeV to 10 TeV.Right panel: Same as left for the γ rays.
S h 2 = 0.3 DM S h 2 = 0.1 DM S h 2 = 0.03 DM S h 2 = 0.01 DM 50.0 52.5 55.0 57.5 60.0 62.5 65.0 67.5 70.0 m S [GeV] FIG.6.95% CL upper limits for the parameter λHS as a function of the DM mass mS.We show the constraints obtained from a reinterpretation of the ATLAS search in Ref.[100] and the projected limits for the HL-LHC and HE-LHC at 14 and 27 TeV, respectively, taken from[73].

FIG. 12 .
FIG.12.Same as in Fig.9but rescaling, in each point of the parameter space, all predicted astrophysical signals (for direct and indirect detection) according to the DM fraction ξ, as it is predicted from the relic density computation.The green shading denotes the region with Ωh 2 > 0.12 (ξ > 1).The orange shaded region shows the region excluded by the direct detection constraints from LZ data.The right panel show a zoom in the region of the resonance, better exposing the configurations compatible with all constraints.
For mS = 31 GeV and mS = 62 GeV, ⟨σv⟩ is dominated by the f f final state; ⟨σv⟩ for mS = FIG. 5. Left panel: Combination of parameters λHS and mS reproducing ΩSh 2 within a fraction (from 1% to 100%) of the measured value of the relic abundance obtained by the Planck experiment, ΩDMh 2 .Here ΩSh 2 is computed solving the full Boltzmann equation taking into account both the elastic and annihilation collision operators (model fBE with QCDB with DRAKE).Right panel: Comparison of the results obtained by solving the Boltzmann equation with both micrOMEGAs and DRAKE [see Eq. (10)].The micrOMEGAs line accounts only for the annihilation processes (nBE); the DRAKE lines assume always fBE for the two extreme QCD phase transition choices, QCDA and QCDB.
[28] panel: Best-fit SED for DM in the SHP model to the GCE found in Ref.[28].We show the data with the statistical and systematic errors (grey band) and the theoretical best-fit (blue curve) with the 3σ uncertainty band.Right panel: Contour plot for the −∆χ 2 defined as −(χ 2 (mS, λHS) − χ 2 min ), where χ 2 min is the chi-square of the best fit.
[62]]panel: Contour regions obtained with a combined fit to cosmic particle flux data (GCE, dSphs and antiprotons).We also show the upper limits on λHS obtained from ATLAS data[100](red dashed line) and the HL-LHC projections for 27 TeV[73](orange dotted line).Direct detection upper limits refer to LZ data[62](brown dashed line) and projections to DARWIN (purple dot-dashed).We also report the region of the parameter space compatible with the observed DM relic abundance via thermal freeze-out (green region) in case of model fBE, including the uncertainty coming from the different choice of the QCD correction approach, labelled with QCD A and QCD B. Right panel: Same as in the left panel but for the region around the resonance mS ≈ m h /2.
[116]4.Same as in Fig.9but with the upper limits derived from CR antiprotons flux data from AMS-02, with different assumptions for the propagation setup.In particular the model labed as 'DIFF' ('INJ') stands for the model DIFF.BREAK (INJ.BREAK) in Ref.[116].Moreover, the cases labeled as 'Unc.' and 'Corr.' refer to the case the uncertainties are taken to be uncorrelated and correlated, respectively.