Reconciling Higgs physics and pseudo-Nambu-Goldstone dark matter in the S2HDM using a genetic algorithm

We investigate a possible realization of pseudo-Nambu-Goldstone (pNG) dark matter in the framework of a singlet-extended 2 Higgs doublet model (S2HDM). pNG dark matter gained attraction due to the fact that direct-detection constraints can be avoided naturally because of the momentum-suppressed scattering cross sections, whereas the relic abundance of dark matter can nevertheless be accounted for via the usual thermal freeze-out mechanism. We confront the S2HDM with a multitude of theoretical and experimental constraints, paying special attention to the theoretical limitations on the scalar potential, such as vacuum stability and perturbativity. In addition, we discuss the complementarity between constraints related to the dark matter sector, on the one hand, and to the Higgs sector, on the other hand. In our numerical discussion we explore the Higgs funnel region with dark matter masses around 60 GeV using a genetic algorithm. We demonstrate that the S2HDM can easily account for the measured relic abundance while being in agreement with all relevant constraints. We also discuss whether the so-called center-of-galaxy excesses can be accommodated, possibly in combination with a Higgs boson at about 96 GeV that can be the origin of the LEP- and the CMS-excess observed at this mass in the bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}-quark and the diphoton final state, respectively.


Introduction
The discovery of a Higgs boson with a mass of approximately 125 GeV at the LHC by both the ATLAS [1] and the CMS [2] collaborations is a milestone for the understanding of the laws of nature. So far, the experimental observation related to the discovered particle, denoted by h 125 in the following, agree with the interpretation of a fundamental scalar particle that behaves according to the prediction of the Standard Model (SM) [3][4][5]. As a result, any model has to incorporate a particle state that resembles a SM-like Higgs boson within the current experimental uncertainties. In contrast to the measurement of the Higgs-boson mass, the measurements of the couplings of h 125 are much less precise, with uncertainties at the level of about ten to twenty percent [4,5]. This leaves room for interpretations of the discovered Higgs boson in models beyond the SM (BSM), in which the theoretical predictions of the properties of h 125 only agree with the SM interpretation at the level of the current experimental uncertainties. In such models, phenomena can be accommodated that cannot be explained in the SM, and the precise measurements of the properties of h 125 will be crucial in order to shed light on the question which of the proposed BSM scenarios could be realized in nature.
Another experimental milestone consists of various indications for the existence of dark matter (DM) through a conjoint of data gathered from, among others, rotation curves of spatial galaxies [6,7], gravitational lensing [8], and the Bullet cluster merger [9]. The Planck -1 -

JHEP10(2021)215
collaboration [10], using the precise map of the cosmic microwave background (CMB), reports the most precise measurement of today's DM relic abundance, (Ωh 2 ) Planck = (0.119 ± 0.003). Hence, the DM sector constitutes about 26% of the energy-matter content of the Universe. Even though there are many indirect indications for the existence of DM via gravitational effects, so far there has not been any direct discovery of a DM particle that could give rise to more information about the properties of DM. The elusive nature of DM has opened up an interesting landscape of BSM theories that can provide one or more DM candidates. One of the most studied scenarios in such SM extensions is the weakly-interacting particle (WIMP), a particle with weak couplings to the SM particles and a mass around the electroweak (EW) scale whose existence could potentially be probed also at present or future colliders. In view of the fact that the DM particle(s) might not be charged under the SM gauge groups, in which case they are also not coupled directly to the quarks and leptons, the possibility of coupling the DM to the SM only via the Higgs sector, often called Higgs portal [11,12], is an interesting scenario. Many extended Higgs sectors provide a (pseudo)scalar DM candidate fitting the WIMP paradigm. However, they are stringently constrained by DM directdetection experiments [13]. A possibility to evade the constraints from direct-detection experiments is given by the fact that the scattering cross sections between the DM and the nuclei can be momentum-suppressed. A particle that naturally has this feature is the so-called pseudo-Nambu-Goldstone boson (pNG) DM [14][15][16][17][18][19][20]. As a result, BSM models that predict the existence of a stable pNG in order to account for the DM relic abundance have recently gained a lot of attention [21][22][23][24][25][26][27][28][29][30].
The most economic way to introduce pNG DM is to extend the SM by a complex singlet field φ S and demanding that the Lagrangian respects a softly broken global U(1) symmetry under which the singlet field is charged [14]. The pNG DM is then given by the imaginary part of φ S , where the global U(1) symmetry prevents the particle from decaying, and the soft U(1)-breaking term gives rise to the mass of the pNG state. The SM has short-comings beyond the fact that it does not provide an explanation for the observation of DM, such that it is also compelling to investigate pNG DM in models with Higgs sectors that compared to the SM contain additional fields together with Φ S . One of such possibilities is that pNG DM can be incorporated into 2 Higgs doublet models (2HDM) (see ref. [31] for a review), in which, in contrast to the pNG DM model with only one Higgs doublet [32], a first-order EW phase transition can be realized [29,33]. Such a transition provides a departure from thermal equilibrium, required to explain the matter-antimatter asymmetry via the mechanism of EW baryogenesis [34,35]. A cosmological first-order phase transition can also source the formation of a stochastic gravitational-wave background that could be observable at future space-based gravitational-wave interferometers, such as LISA [29]. Moreover, we emphasize that in the S2HDM the presence of the second Higgs doublet gives rise to two additional neutral and two charged scalars. For DM masses comparable to the masses of these additional states, new annihilation processes can become important for the prediction of the DM relic abundance. As a result, the predicted relic abundance for a certain DM mass can differ substantially between the S2HDM and the simpler model with only one Higgs doublet, and new parameter regions can become physically viable, where the corresponding parameter space of the model with only one Higgs double predicts a too large relic abundance [36].

JHEP10(2021)215
In addition to the above mentioned phenomenological reasons to consider a model with a second Higgs doublet field, one should also note that supersymmetric extensions of the SM, in which the hierarchy problem can be addressed [37][38][39], require the existence of at least two Higgs doublet fields φ 1,2 in order to account for the masses of all quarks and leptons. Also the most commonly studied solutions to the strong CP-problem incorporating the so-called QCD axion require the presence of two doublet fields [40]. Moreover, new axially coupled U (1) interactions, resulting in extra gauge bosons weakly coupled to standard model particles and which behaves very much as an axion-like particle, provide a possible bridge to a new dark sector and also demand an additional electroweak doublet [41]. Other models to solve the hierarchy problem rely on a unification of the gauge interactions and the fact that h 125 arises as a (composite) pNG, but where also additional (potentially stable) pNG can be present [42]. Such models could resemble at low energy a model with two Higgs doublets [43][44][45][46].
In this paper we study a singlet extension of the 2HDM (S2HDM), where the real part of the singlet field φ S gives rise to a third Higgs boson, and the imaginary part of φ S gives rise to the pNG DM, as discussed above, while the terms of the scalar potential incorporating the doublet fields φ 1,2 are identical to the 2HDM with softly broken Z 2 symmetry. In total, the physical scalar particle spectrum consists of three CP-even Higgs bosons h 1,2,3 that are mixed with each other, two CP-odd states A and χ that do not mix and where χ is the stable DM candidate, and finally the charged Higgs bosons H ± . We focus on the S2HDM with type II Yukawa structure and the parameter space that gives rise to the so called Higgs funnel scenario, i.e. resonant DM annihilation via s-channel diagrams mediated by the SM-like Higgs boson h 125 . The possibility of relatively light pNG DM in a doublet extension of the SM suggest an interesting interplay of collider phenomenology and astrophysics that can be constrained by various experimental requirements coming from flavour physics, electroweak precision observables, searches for additional scalars and measurements of the properties of h 125 . Additionally, there are experimental constraints related to the presence of the DM candidate. In particular, the limitation imposed by the fact that a too large relic abundance after thermal freeze-out would overclose the universe and indirect-detection limits coming from the observation of dwarf spheroidal galaxies (dSph) by the Fermi-LAT space telescope [47] play an important role. We also take into account several theoretical constraints that must be imposed in order to ensure the validity of the perturbative treatment of the theory and the stability of the EW vacuum. While the S2HDM was already studied in regards to the DM phenomenology and its cosmological history in refs. [29,36,46], a careful treatment of the S2HDM taking into account the large number of constraints has not been carried out yet. In our analysis we will demonstrate that the combined consideration of the experimental and theoretical constraints is crucial in order to make reliable predictions for the phenomenology of the S2HDM.
In this context, we explore benchmark scenarios featuring pNG DM in the mass range 40 GeV ≤ m χ ≤ 80 GeV. This region is particularly relevant from the experimental point of view, since it belongs to parameter space of the S2HDM where it is possible to accommodate a sizable fraction of the DM relic abundance, and, whenever this is the case, the presence of the DM candidate is currently, and even more so in the new future, probed by DM -3 -

JHEP10(2021)215
indirect detection experiments. As a result, important limitations on the parameter space of the S2HDM will arise. In this regard, it is interesting to note that the corresponding parameter space is also suitable to realize the excess of gamma rays from the galactic center observed by the Fermi Large Area Telescope (LAT) [48,49]. It has been argued that these observations could be originated by DM annihilations in the galactic center [50][51][52][53][54][55][56] where a large concentration of DM is expected to reside [57,58]. However, they are also consistent with an unresolved population of thousands of millisecond pulsars in the Galactic bulge [59][60][61]. At the same time, the Alpha Magnetic Spectrometer (AMS) [62], onboard the International Space Station, reported an excess over the expected flux of cosmic ray antiprotons consistent with DM annihilating into b-quark pairs with a similar range of DM masses [63][64][65][66][67][68]. We will address the question whether the DM candidate of the S2HDM can account for the two cosmic-ray excesses, potentially in combination with a Higgs boson at roughly 96 GeV that could give rise to the so-called LEP excess in the bb final state and the CMS excess in the diphoton final state, which was already investigated in a singlet-extension of the SM featuring pNG DM in ref. [25]. Therein, it was found that the CMS excess requires new charged particles in order to account for a sufficiently strong signal. In contrast, here we will follow the results of ref. [69] obtained in the Next-to 2HDM (N2HDM) and in models with supersymmetry [70][71][72][73][74][75]. It was shown that the presence of a second Higgs doublet in addition to a singlet scalar field allows for an explanation of both collider excesses without having to rely on new charged states that appear in the loops of the loop-induced coupling of the possible Higgs boson at 96 GeV in order to enhance its diphoton rate. This paper is organised as follows. In section 2, the model is presented. In section 3, we describe the relevant experimental and theoretical constraints that we apply to its parameter space. In section 4, we describe the genetic algorithm that was used to scan the parameter space and to determine the parameter points that pass the various theoretical and experimental requirements. In section 4.1, we explore the Higgs funnel region after imposing the previously described constraints and disregarding the explanation of the excesses at LEP and CMS, whereas in section 4.2 we additionally demand that the collider excesses are accommodated. Finally, we conclude in section 5.

The S2HDM
The scalar sector of the S2HDM consists of two SU(2) doublets and a complex gauge singlet field, which can be expressed as where the pseudoscalar component χ gives rise to the DM candidate of the model. Assuming the absence of explicit CP violation, the scalar potential of the S2HDM is given by

JHEP10(2021)215
Here the terms that exclusively involve the doublet fields are identical to the scalar potential of the 2HDM, where a Z 2 symmetry defined by the transformations φ 1 → φ 1 , φ 2 → −φ 2 and φ S → φ S is only softly broken by the terms proportional to µ 12 . After the generalization of the Z 2 to the Yukawa sector, which is unchanged in the S2HDM compared to the 2HDM, the Z 2 symmetry gives rise to the absence of flavour-changing neutral currents at tree level. Depending on the assigned charges of the fermions, this results in the usual four Yukawa types (see ref. [31] for details). We will focus on the so-called type II, familiar from supersymmtric BSM scenarios, in which φ 2 is only coupled to up-type quarks, and φ 1 is only coupled to down-type quarks and the charged leptons. The remaining terms of the scalar potential V involve the singlet field φ S and respect a global U(1) symmetry, except for the term proportional to µ 2 χ . This term softly breaks the U(1) symmetry, thus providing a non-zero mass for the pNG DM.
Without loss of generality, the field configuration of the vacuum can be expressed as where we made use of the fact that redundant degrees of freedom related to the gauge symmetries can be removed via the gauge transformations. We will focus on the case in which the EW symmetry is broken by non-zero values of v 1 and v 2 , and an accidental Z 2 symmetry, under which ρ S changes the sign, is broken by v S > 0. The charge-breaking vev v C , the CP-breaking vev v CP , and v DM are considered to be vanishing, noting that a non-zero value of v DM would give rise to decays of the DM candidate χ. In our numerical analysis, we verify for each parameter point whether there exists a global minimum of the potential with v 1 , v 2 , v S > 0 and v C , v CP , v DM = 0. Otherwise, we remove such a point, since it potentially features an EW vacuum that is short-lived compared to the age of the universe (see section 3.1 for details). In order to make a connection to the SM and the 2HDM we define the parameters v 2 = v 2 1 + v 2 2 ∼ (246 GeV) 2 and tan β = v 2 /v 1 . Assuming the EW vacuum configuration as described above, the CP-even fields ρ 1,2,S mix, giving rise to the mass eigenstates h 1,2,3 , where throughout this paper the mass hierarchy m h 1 < m h 2 < m h 3 will be assumed. The mixing in the CP-even sector can be written in terms of an orthogonal transformation given by a matrix R, such that where −π/2 ≤ α 1 , α 2 , α 3 ≤ π/2 are the three mixing angles, and we use the short-hand notation s x = sin x, c x = cos x. The charged scalar sector remains unchanged compared to the 2HDM. It contains two physical charged Higgs bosons H ± with mass m H ± and the charged Goldstone bosons related to the gauge symmetries. The pseudoscalar components σ 1 , σ 2 and χ form a neutral Goldstone boson and two physical states A and χ with masses m A and m χ , respectively. Here it is important to note that the remnant Z 2 symmetry that is present when v DM = 0, preventing the DM candidate χ from decaying, also forbids the -5 -

JHEP10(2021)215
mixing between χ and A. Hence, the pseudoscalar A has effectively the same couplings to the fermions as the one of the 2HDM.
Given the definitions of the parameters as defined above, it is possible to replace most of the parameters of the scalar potential shown in eq. (2.2) by more physically meaningful parameters. In our numerical analysis, we will sample the parameter space of the type II S2HDM in terms of the parameters The relations between the parameters shown in eq. (2.5) and the Lagrangian parameters of the potential in eq. (2.2) are given in appendix A.

Constraints
In this section we briefly discuss the constraints on the parameter space of the S2HDM that we applied in our analyses. In each case, we also illustrate the impact of the constraint in order to give an impression on their relevance for our numerical discussion. Some of the constraints, such as the ones arising from demanding a stable EW vacuum, are similar to the corresponding constraints known from the Next-to 2HDM (N2HDM) [76]. However, there are also important differences which we will point out below.

Theoretical constraints
The parameters that appear in the scalar potential of the S2HDM are subject to important theoretical constraints. These constraints ensure the stability of the EW vacuum for a given parameter point, and they exclude parameter values which give rise to parameter points that could not be treated perturbatively.
Boundedness-from-below. We apply bounded-from-below (BfB) conditions on the tree-level scalar potential, which determine whether the potential is bounded from below along all field directions. Due to the fact that the quartic part of the potential V is unchanged compared to the N2HDM, we can apply the same conditions that were found for the N2HDM [76,77]. We exclude all parameter points from our analyses that do not feature a scalar potential that is BfB. It was shown in ref. [78] that large loop corrections can transform a bounded tree-level 2HDM potential into an unbounded one, potentially destabilizing the EW vacuum. This effect is expected to be also present in the S2HDM, such that our tree-level analysis of the boundedness could permit potentially unphysical parameter points. However, the possibility of loop corrections changing the boundedness of the potential was shown to be present only in regions of the parameter space with splittings between m A , m H and m H ± larger than ∼ 250 GeV, where consequently large quartic couplings are present [78], which then give rise to the corrections. In our analysis, we demand an upper limit of 200 GeV on the splitting of the heavy Higgs-boson masses compared to the mass scale M defined in eq. (2.5) (see also the discussion below), such that we expect that the boundedness of the potential, and therefore the stability of the EW vacuum, are not to be severely affected by the loop corrections. EW vacuum stability. In the next step, we demanded that the EW minimum as described in section 2 is the global minimum of V , such that no vacuum decay into other unphysical minima is possible, and the stability of the EW vacuum is guaranteed. To test whether the EW minimum is the global one, we first determined all extrema of V by solving the stationary conditions where we used the code Hom4PS-2 [79] to solve the system of polynomial equations. For each extrema we calculated the value of V in this point of field space. One can conclude that, if for any of the extrema the value of V is smaller than for the field values of the EW vacuum, the EW minimum is not the global minimum of V . In this case, the EW vacuum is potentially short-lived compared to the age of the universe, such that the corresponding parameter point might be unphysical, and we rejected it from the analyses. 1 Perturbative unitarity. In order to verify whether a perturbative treatment of the model is valid for a given parameter point, we applied the so-called tree-level perturbative unitarity constraints. We derive these constraints by calculating the scalar 2 × 2 scattering matrix in the high-energy limit, in which only the quartic contact interactions are relevant. The precise form of the conditions is given in appendix B for completeness. They set upper limits on the absolute values of the parameters λ i and combinations thereof, such that they are particularly relevant when there are large mass splittings between the heavy doublet states A, H ± and one of the scalars h i . Due to the fact that compared to the N2HDM the only additional degree of freedom is the CP-odd component of the singlet field φ S , the perturbativity conditions are in most parts very similar to the N2HDM conditions [76]. However, an important difference is that an additional condition on the singlet self-coupling of the form |λ 6 | < 8π appears. In addition, the constraints related to scattering amplitudes involving the singlet field components and the field components of the doublet fields (see eq. (B.14)) are modified with respect to the N2HDM.
Energy scale dependence of the theoretical constraints. Both the perturbative unitarity constraints and the BfB conditions are in many analyses of the 2HDM or its extensions applied exclusively at a certain energy scale. However, it is known that the model parameters obtain an intrinsic energy scale dependence due to radiative corrections, which is governed by their evolution under the running group equations (RGE). It is therefore possible that even though at the initial scale, here assumed to be µ = v such that the input scale of the parameters corresponds to the EW scale, a parameter point passes the theoretical constraints, the point becomes unphysical at larger energy scales µ > v. Due to the fact that the pertubativity conditions allow for values of |λ i | > 1, the 1 A (zero temperature) calculation of the lifetime of an unstable EW vacuum shows that in some cases the EW vacuum can be considered to be sufficiently long-lived, even though there are deeper minima present, such that a parameter point with a non-global EW minimum could still be viable (see ref. [80] for an N2HDM analysis). However, in such cases it is still unclear whether the universe would have adopted the (meta-stable) EW vacuum at some point within the thermal history of the universe, or would have rather transitioned into a deeper unphysical minimum. The analysis of the thermal history of the scalar potential of the S2HDM is beyond the scope of this paper (see ref. [33] for an N2HDM analysis), such that we demand the most conservative constraint, i.e. excluding all parameter points for which the EW minimum is not the global minimum of the potential. energy range v ≤ µ ≤ µ v in which the theoretical constraints are fulfilled can be small, with values of µ v even smaller than the energy scales that are probed at the LHC. In consequence, we apply the previously described theoretical constraints taking into account the energy-scale dependence of the parameters, utilizing the two-loop β-functions of the S2HDM and demanding that the theoretical constraints are respected up to a certain energy scale µ v . The β-functions for the S2HDM were obtained with the help of the public code SARAH [81,82], solving the general expressions known in the literature [83][84][85]. We also calculated the β-functions with the code PyR@TE 3 [86] to be able to cross check the expressions and found exact agreement. We discarded a parameter point when the scale µ v at which the scalar potential becomes unbounded or at which the perturbative unitarity constraints are violated is below 1 TeV, which was also chosen as the upper limit on the Higgs-boson masses in the numerical discussion (see section 4).

Experimental constraints
The S2HDM offers a rich phenomenology that can be probed experimentally by various means. The corresponding experimental (null)-results give rise to numerous constraints that have to be taken into account. We start by discussing the constraints related to the Higgs sector of the model. Subsequently, we describe the manner in which the constraints from measurements from DM experiments were taken into account.

Searches for additional scalars and properties of h 125 .
Regarding the Higgs phenomenology of the model, we used the public code HiggsBounds v.5.9.0 [87][88][89][90][91][92] to test the parameter points against a large number of cross-sections limits from direct searches for Higgs bosons at LEP, the Tevatron and the LHC. For each Higgs boson, HiggsBounds selects the potentially most sensitive experimental search based on the expected limits. For the selected searches, the code then compares the predicted cross sections against the observed upper limits on the 95% confidence level and excludes a parameter point whenever the theoretical prediction lies above the experimental limit for one of the Higgs bosons.
Regarding the discovered Higgs boson at 125 GeV, we use the public code HiggsSignals v.2.6.1 [93][94][95][96] to verify whether an S2HDM parameter point features a particle h i that resembles the properties of the discovered particle h 125 within the experimental uncertainties. HiggsSignals performs a χ 2 -analysis confronting the predicted signal rates against the experimentally measured signal rates. In our more general parameter scan discussed in section 4.1, we applied as constraint that the resulting χ 2 value (called χ 2 125 in the following) fulfills χ 2 125 ≤ χ 2 SM,125 + 5.99, where χ 2 SM,125 = 84.41 is the fit result assuming a SM Higgs boson at 125 GeV, and where the allowed penalty of 5.99 corresponds to a 95% confidence interval for two-dimensional parameter distributions. 2 In section 4.2, in which we aim for accommodating the collider excesses observed at about 96 GeV, we combine the value of χ 2 125 obtained from HiggsSignals with a value χ 2 96 that quantifies the fit to the excesses. The precise criterion applied will be given in section 4.2.
HiggsBounds and HiggsSignals require as input effective coupling coefficients, which are defined as the couplings of the physical scalars normalized to the coupling of a SM Higgs 2 See ref. [96] for details on the interpretations of the χ 2 analysis of HiggsSignals. boson of the same mass. With the help of these coupling coefficients, the codes compute the relevant cross sections for the scalars by rescaling the predictions for a hypothetical SM Higgs boson. In the S2HDM, the coupling coefficients can be expressed in terms of tan β and (for the states h i ) in terms of the mixing angles α i . The precise expressions are identical to the N2HDM expressions and can be found in ref. [76]. Moreover, the user has to provide the branching ratios of the Higgs bosons. We calculated these in two steps. First, we used the public Fortran code N2HDECAY [76,[97][98][99][100] implemented in the anyhdecay C++ library to calculate the decay widths of h i , A and H ± for decays into SM particles and for cascade decays with one or two Higgs bosons in the final state. 3 In a second step we calculated the decay widths for the invisible decay into a pair of χ as described below (see eq. (3.1)). We finally divided each partial decay widths by the total widths to obtain the branching ratios for each possible decay mode.
In addition to the global constraints on the measured signal rates of h 125 , the S2HDM can also be probed via possible decays of h 125 into a pair of DM particles χ with a mass of m χ < 125/2 GeV. At leading order, the partial decay width of the invisible decay is given by The most recent upper limit on the branching ratio of the invisible decay BR inv of h 125 has recently been reported by ATLAS and is given by BR inv < 0.11 at the 95% confidence level [101]. We applied this limit as additional constraint complementary to the HiggsSignals analysis. However, as will be demonstrated in section 4.1, in most cases parameter points with sizable values of the corresponding branching ratios BR inv are already excluded by the global constraints on the measured signal rates of h 125 , since the additional decay mode h 125 → χχ suppresses the ordinary decays of h 125 into SM final states.
Electroweak precision observables. Further constraints originating from the presence of the BSM Higgs bosons that have to be taken into account are the ones related to EW precision observables (EWPO). Since the S2HDM extends the SM particle content exclusively by scalar states, one can to a very good approximation apply the formalism of the oblique parameters S, T and U [102,103]. These parameters are theoretically defined by loop corrections to the gauge-boson self-energies, in which the BSM particles appear in the loops, giving rise to modifications of the values of the oblique parameters compared to the SM. In order to predict the oblique parameters, we applied the general expressions at the one-loop level from refs. [104,105] to the S2HDM. Experimentally, S, T and U are constrained via global fits to the EWPO, where we utilize here the results (including their uncertainties) found in ref. [106]. In 2HDM-like extensions of the SM, the most sensitive parameter is the T parameter, whereas the modifications of the U parameter in practically all cases are orders of magnitude smaller than the experimental sensitivity, and we explicitly checked this to hold in the S2HDM. 4 We therefore performed a two-dimensional χ 2 test regarding S and T , written as χ 2 ST in the following, and discarded parameter points for which the predicted values were not in agreement with the experimental fit result [106] at the 95% confidence level. This gives rise to the requirement χ 2 ST ≤ 5.99. The T parameter is sensitive to the breaking of the custodial symmetry. As a result, one finds strong exclusions when there is a sizable mass splitting between the states A, H ± and (depending on the doublet-admixture) one of the heavy CP-even state h 2 or h 3 .
Flavour-physics observables. Also the theoretical predictions for flavour-physics observables are modified compared to the SM via contributions from the additional Higgs bosons of the S2HDM. In particular, the presence of the charged Higgs bosons gives rise to robust constraints in the parameter plane of tan β and m H ± . Since there are no public results for the theoretical predictions in singlet extensions of the 2HDM for some of the most relevant flavour observables, we simply applied hard cuts on the ranges of tan β and m H ± in our numerical analysis, where the cuts were determined by assuming that the exclusion regions known from the 2HDM are not severely modified by the presence of the additional field of the S2HDM, which we expect to be the case due to the singlet nature of this field. Consequently, working in the type II S2HDM, we set lower limits of tan β > 1.5 and of m H ± > 600 GeV in order to not be in conflict with constraints from radiative and (semi-)leptonic B meson decays and from their mixing frequencies [106]. 5 Dark matter observables. We now turn to the experimental constraints that are related to the presence of the dark matter candidate χ. The most important limitation arises from the fact that a too large relic abundance of χ after thermal freeze-out would overclose the universe. The currently most precise measurement of today's DM relic abundance Ωh 2 is given by surveying the cosmic microwave background by the Planck satellite, leading to a measurement of (Ωh 2 ) Planck = (0.119 ± 0.003) [10]. We will use this value as an upper limit on the relic abundance of χ in our analysis, taking into consideration that in case the relic abundance of χ is smaller than (Ωh 2 ) Planck there is room for additional (particle or astrophysical) contributions to the relic abundance. We focus the analysis on the Higgs funnel region with DM masses of 40 ≤ m χ ≤ 80, where there are good prospects to be able to explain most (or all) of the observed DM relic abundance via the thermal freeze-out of χ [14,23,25,29,36]. For the theoretical prediction of the relic abundance, we wrote an S2HDM modelfile for the Mathematica package FeynRules v.2 [109][110][111], which we utilized to obtain a CalcHEP [112] input for the public code MicrOMEGAs v.5 [113] written in C and Fortran. With this input, MicrOMEGAs is capable of calculating the relic abundance and the freeze-out 4 We found that at the one-loop level the theoretical predictions for S, T and U in the S2HDM and the N2HDM (given the same values of m h i , mA and m H ± ) are identical, because they do not depend on the additional state χ of the S2HDM as long as v DM = 0. 5 A more recent result suggests a lower limit of m H ± > 800 GeV in the type II 2HDM from the measurement of the radiative B meson decay [107], whereas ref. [108] claims that theoretical uncertainties might have been underestimated in the literature, potentially giving rise to a weaker lower limit. We emphasize that the conclusions drawn from our numerical analysis do not depend on the precise value of the lower limit chosen for m H ± .

JHEP10(2021)215
temperature, where for the computation of the annihilation cross section all 2 × 2 processes and also processes with off-shell vector bosons in the final state are taken into account. As already pointed out in section 1, one of the attractive features of the S2HDM is that due to the pNG nature of the DM particle the cross sections for the scattering of χ on nuclei vanish at leading order in the limit of vanishing momentum transfer [36], such that at this order direct-detection experiments are not sensitive to the presence of χ. In addition, it was shown in models with a single Higgs doublet field and a complex singlet field that the loop contributions to the direct-detection cross sections are small, and the predicted direct-detection scattering cross sections remain far below the current (and near future) sensitivity of direct-detection experiments [21,22,28]. In the type II S2HDM the masses of the additional doublet particle states H(= h 2 or h 3 ), A and H ± are required to be substantially larger than the DM masses m χ considered in our analysis (see discussion above), such that we can safely assume that the relevant loop corrections to the DM-nuclei scattering cross sections are captured by the pNG DM model with only one Higgs doublet. One should note also that for light DM an additional suppression of the scattering cross section given by the factor (m χ /m h i ) 4 is present, which reflects the pNG nature of χ and the fact that the loop corrections vanish in the limit m χ → 0 [17]. Consequently, in our scenario the additional loop corrections arising from the presence of the second Higgs doublet are even smaller than the corrections known from the pNG DM model with one Higgs doublet, and thus there are no relevant constraints from direct-detection experiments that have to be taken into account in our analysis (see also discussions in refs. [29,36]).
On the other hand, constraints from DM indirect-detection experiments are important, in particular in the Higgs funnel region investigated here, in which χ mainly annihilates into b quark pairs, typically via h 125 in the s-channel. The most stringent constraints on the annihilation cross sections of DM come from the observation of dwarf spheroidal galaxies (dSph) by the Fermi-LAT space telescope [47]. In order to account for these constraints, we used FeynRules to generate UFO [114] model files for the S2HDM, which were then used as input for the public code MadDM v.3 [115,116]. MadDM is a plugin for MadGraph5_atMC v.3.1.1 [117] that can be used to compute the relevant velocity-averaged annihilation cross sections σv rel bb , and to subsequently compare the theoretical predictions to the upper limits on the velocity weighted cross section for DM particles annihilating into bb final states from the Fermi measurements of gamma rays from dSph at the 95 % CL. 6 The Fermi-LAT collaboration utilizes a likelihood analysis to fit the spectral and spatial features of dSphs to obtain upper limits on the annihilation cross section as a function of the DM mass [47]. The analysis accounts for point-like sources from the latest LAT source catalog, models the galactic and isotropic diffuse emission, and incorporates uncertainties in the determination of astrophysical J-factors, which depend on both the DM density profile and the distance. The observed limits are sensitive to the determination method of the J-factors. In ref. [47] an evaluation of the uncertainties arising from targets lacking measured J-factors was performed. Using only predicted J-factors for the whole sample weakened the observed limits by a factor of about 2 to 3, depending on the choice of J-factor uncertainty, with respect to the limits obtained by using both predicted and measured J-factors. Considering these uncertainties will be important for the discussion of the tension between the constraints coming from dSph and the gamma-ray excesses and anti-protons measured from the galactic center, as will be demonstrated in section 4.1.
For the comparison between the predicted annihilation cross section and the Fermi bounds from dSph observations, we rescaled (when not explicitly said otherwise) the cross sections with a factor in order to account for the suppression of today's annihilation cross section of χ due to the smaller number density when the relic abundance of DM is not made up completely out of χ. 7 We also point out that the velocity-averaged annihilation cross sections can be considered here to be velocity-independent in the non-relativistic limit to a very good approximation, since in our scan range of m χ they are dominantly generated via diagrams with s-channel exchange of either h 1 or h 2 [118]. Nevertheless, we calculated σv rel with different relative velocities v rel for the comparison against the Fermi-LAT dSph constraints, on the one hand, and for the comparison against the preferred regions regarding the gammaray and the anti-proton excesses, on the other hand. In both cases we used the default values of MadDM, which are v rel = 2 · 10 −5 for the DM in dSph and v rel = 10 −3 for DM in the center-of-galaxy as relevant for the excesses. In agreement with our expectation, the differences of the annihilation cross sections for the two values of v rel stayed below a few percent and are not relevant for our discussion.

Numerical analysis
As was already discussed in section 1, we divide our numerical analysis of the type II S2HDM into two parts. In the first part discussed in section 4.1, we will demonstrate in a broad parameter scan how the Higgs funnel region with 40 GeV ≤ m χ ≤ 80 GeV is affected by the various theoretical and experimental constraints. Here the DM particle χ is the lightest BSM state, and h 125 = h 1 is the lightest of the three CP-even Higgs bosons h i . We will describe in detail the predictions for the DM relic density and its interplay with the Higgs sector of the model. In addition, we investigate whether the annihilation of χ in this scenario could give rise to the cosmic-rays anomalies from observations of the spectra of cosmic rays coming from the center of the galaxy. We emphasize at this point that due to the large mass gap between the DM mass m χ studied here and the masses of the heavy scalar states h 3 , A and H ± in the type II S2HDM, the predictions for the DM relic abundance and today's DM annihilation cross section mainly depend on the couplings

JHEP10(2021)215
of χ to the SM-like Higgs boson and (when present) the light singlet-like Higgs boson. Accordingly, the properties of the DM sector will be similar compared to the predictions from the pNG DM model with only one Higgs doublet, because additional annihilation processes involving the heavier states (see also the discussion in section 1) do not play a role. However, differences between both models can still arise due to the richer mixing patterns of the states h i in the S2HDM, where the mixing angles α 1,2,3 enter the couplings of h i to χ.
In the second part of our analysis, discussed in section 4.2, we focus on the parameter space in which at the same time the collider excesses at about 96 GeV could be accommodated. Consequently, here the presence of a singlet-like Higgs boson h 96 = h 1 with m h 1 = 96 GeV is enforced as an additional constraint on the parameter space. As a result, the SM-like Higgs boson h 125 is the second lightest Higgs boson h 2 , and its mixing with h 96 is subject to the constraints from the LHC measurements of the signal rates of h 125 . Going beyond the discussion of the collider phenomenology and the excesses at 96 GeV, we will illustrate in detail how the presence of h 96 has also important consequences for the DM phenomenology in the Higgs funnel, in particular giving rise to a second s-channel contribution to the thermal freeze-out cross section and today's annihilation cross section relevant for DM indirect-detection experiments.
In both parameter scan presented in the following, we sampled the multi-dimensional parameter space of the model utilizing a genetic algorithm. In contrast to random or uniform (grid)-scans of the model parameters, a genetic algorithm has the advantage that it focuses on the relevant parameter region by minimizing a so-called loss function, which has to be suitably defined in each case. The definition of the loss functions used in both parts of our analysis will be given in section 4.1 and section 4.2. Apart from the loss function, the properties of the genetic algorithm applied were in large parts identical in both scans. For the interested reader we briefly describe the main design choices here, where we made use of the public python package DEAP [119] to perform the algorithm.
The algorithm starts by generating an initial sample (also called population) of 50 000 parameter points. Each parameter point (also called individual) is defined by a list of 14 numbers (also called attributes or genes), where each number of this list defines a value of one of the model parameters within a given parameter range. The population is then subject to an evolution including the three steps: selection, mating and mutation. These three steps are performed in a loop for a total number of N cycles (also called generations), such that each cycle gives rise to a new population of parameter points with (desirably) better fitnesses. The fitness of each individual is defined by the corresponding value of the loss function: the smaller the value of the loss function given the parameter values of an individual, the better is the fitness of the individual.
The first step of each cycle, i.e. selection, determines which of the individuals of the population are allowed to take part in the following two steps, i.e. mating and mutation. As a selection function we used the so-called tournament selection with size three. This function selects the individual with the best fitness from three randomly picked individuals of the population. In total 50 000 individuals are selected in this way (where each individual was allowed to be selected more than once) and these then proceed to the mating stage.

JHEP10(2021)215
Since the selection is based on the fitness values, individuals with better fitness have a higher chance of producing new individuals (called offspring).
For the mating process, we divided the selected individuals into two distinct groups, and then we performed a uniform crossover of pairs of individuals from each group. A uniform crossover creates two child individuals from each pair of parent individuals, where the child individuals are defined by swapping the attributes of the two parent individuals, in our case according to a probability of 0.2. Hence, the two parent individuals produce two offspring individuals which have on average 20% of the attributes from one parent and 80% of the attributes from the other parent. In addition, we included a so-called mating probability of 0.8, such that for 20% of the pairs of parent individuals no mating was performed and the parent individuals were just kept in the population without changing their attributes.
Afterwards, the mutation stage is performed, which modifies some of the individuals of the offspring via a randomized function, potentially giving rise to new individuals with good fitness values that belong to so far unexplored parameter regions. As a mutation function we applied the so-called float uniform mutator function with a mutation probability of 0.2. This function multiplies the attributes of an individual with a random number between 0.8 and 1.2 according to a probability of 0.1. As a result, 20% of the individuals of the offspring are mutated, and the mutations modify on average 10% of the attributes of such individual.
At the end of each cycle, we replace the initial population with the offspring and enter a new cycle, until either an individual is found that corresponds to a value of the loss function below a certain threshold, or until the maximum number of cycles is reached. Since it is possible that the individual in the parent population with the best fitness would be lost when the population is replaced, we append this best-fit individual to the offspring population in order to ensure that it always survives the complete cycle. Finally, when the algorithm has completed, we save the parameter point with the best fitness. Accordingly, the above described algorithm is performed as many times as the number of desired parameter points in the final sample.
For the two scans discussed in section 4, we compared the performance of the genetic algorithm to the one of a random scan over the free parameters using a flat prior. For a machine-independent estimate of the performances of both algorithms, we chose the number of evaluations of the loss function L (see eq. (4.2)) that is required until a parameter point featuring a value of L below a certain threshold is found. We found for the first scan discussed in section 4.1 that, on average, the genetic algorithm succeeds in finding a value of L < 90 with roughly 60% to 70% fewer evaluations of L compared to the random scan, such that the improvment is only moderate. For the second scan discussed in section 4.2, in which L receives an additional term, our computations indicate that the genetic algorithms outperforms the random scan drastically. Here we found that using the genetic algorithm the average number of evaluations of L in order to find a parameter point with L < 150 was approximately 35 times smaller than using a random scan. Since in this scan the parameter points with the desired features with regards to the collider excesses (see section 4.2 for details) require values of L that are even smaller than L = 150, we conclude that the usage of the genetic algorithm was a vital piece of our numerical analysis. The reason for the fact that the genetic algorithm performs so much better in the second scan, whereas the -14 -

JHEP10(2021)215
improvement was only moderate in the first scan, can be attributed to the fact that the simultaneous minimization of the values of χ 125 (see section 3.2) and the value χ 96 (defined in section 4.2), which quantifies the fit to the collider excesses, requires additional relations between the mixing angles α i and tan β, which the genetic algorithm is able to find more quickly by successively adjusting the parameters of the points with the lowest values of L that have been found in the previous generation.

pNG DM in the Higgs funnel region
In order to explore the Higgs funnel region, we scanned the parameter space of the S2HDM within the parameter ranges Thus, this condition on ∆M max ensures that the masses of the heavy doublet-like states A, H ± and H = h 2 or = h 3 are not further than 200 GeV away from the mass scale M . As explained in section 3.1, and as will also be demonstrated in the following, this condition excludes parameter points that have a very small energy range v ≤ µ ≤ µ v in which the parameter points fulfill the theoretical constraints, with potentially µ v 1 TeV. The lower limits of tan β ≥ 1.5 and m H ± ≥ 600 GeV exclude parameter points that are potentially in conflict with constraints from flavour-phyiscs observables. The mass hierarchy of the CP-even Higgs bosons h i is fixed such that h 125 = h 1 is the lightest one. Their mixing angles α i are scanned over the theoretically possible range, where it should be noted that their values are strongly constrained by the measurements of the signal rates of h 125 , as will also be demonstrated below. The vev of the singlet field v S is allowed to take on values up to 1 TeV, which coincides with the upper value chosen for the masses of the heavier BSM states H ± , A and h 2,3 . If we would have allowed for larger values of v S and M , the heavy states could acquire also larger masses and decouple from the lighter states h 1 and χ. However, we wanted to focus on the parameter space region in which the collider constraints from direct searches at the LHC play a role, such that we limited our scan to the case in which all particle states could be produced (and discovered) at the LHC.
The scan points that we will present were obtained in a two step procedure. In the first step we applied the genetic algorithm as described before in order to find parameter points that minimize the loss function

JHEP10(2021)215
divided by the experimentally observed upper limit (see section 3.2 for details). As a result, parameter points featuring a value of r HB obs > 1 should be rejected, and the second term in the loss function quantifies the penalty of this requirement. The factor 100 is included in order to enhance the importance of this exclusion in terms of the loss function compared to the values of χ 2 125 , thus making sure that all parameter points with low values of the loss function have r HB obs < 1 and are consequently not excluded by direct searches. Finally, the third term is a huge constant C that is added when a parameter point does not fulfill the theoretical constraints at the initial energy scale µ = v, or when the constraints from the EWPO are not fulfilled. With this definition of the loss function, the genetic algorithm finds parameter points that pass the theoretical constraints, the constraints from the collider experiments and the EWPO.
In a second step, all the parameter points found with the genetic algorithm were subject to the remaining constraints: according to the discussion in section 3.1, we applied the theoretical constraints for scales µ > v and verified whether they are fulfilled up to at least µ = 1 TeV. In addition, we verified that, regarding the SM-like Higgs boson, we have ∆χ 2 125 = χ 2 125 − χ 2 SM,125 ≤5.99 and BR inv < 0.11, and, regarding the DM candidate, that the predicted relic abundance is not larger than the Planck value, i.e. Ωh 2 ≤ (Ωh 2 ) Planck . We also ensured that the constraints from the indirect-detection experiments from the observation of dSph are respected. The DM observables were not taken into account already in the definition of the loss function, because the computation of the relevant theoretical predictions were the most time-consuming part of the analysis, such that it was much more efficient to perform these computations only for the parameter points that otherwise passed all the other theoretical and experimental constraints.
As was already mentioned before, the main purpose of this analysis is to illustrate the combined impact of the various constraints on the model parameters. In particular, we will point out which of the constraints give rise to limitations on which subset of parameters, and whether the constraints cover similar or clearly distinct regions of the S2HDM parameters. In the following, we start the discussion with the theoretical constraints that were applied according to the discussion in section 3.1. In the next step, we examine the impact of the collider constraints by taking into account both the constraints from direct searches and from the constraints on the properties of h 125 (see section 3.2). Finally, we consider the physics related to the DM candidate χ, and how its properties are interconnected to the Higgs sector.
In order to analyze the impact of the theoretical constraints, we show in figure 1 the parameter points with the colour coding indicating the energy scale µ v until which the theoretical constraints are respected. We remind the reader that all parameter points fulfill the theoretical constraints at the initial scale µ = v. All points for which µ v < 1 TeV are shown in grey. We performed the RGE running up to µ = 100 TeV, such that points that have µ v = 100 TeV (yellow points) are potentially valid up to much higher energy scales. We find that the relevant constraint that give rise to the low values of µ v are in most cases the tree-level perturbative unitarity constraints. These constraints effectively provide upper limits on the absolute values of the quartic scalar couplings λ i and combinations thereof (see also appendix B). It is therefore easy to understand why they are more severe in region of -17 -

JHEP10(2021)215
parameter space with relatively large splittings between the masses of the heavy BSM states and the mass scale M , since such splittings are induced by large absolute values of λ 1,2,3,4,5 (see also appendix A). Moreover, for obvious reasons also the energy scale dependence of the quartic couplings is stronger when their absolute values are larger. As a result, points with large mass splittings, which potentially were already on the edge of being excluded via the tree-level perturbative unitarity constraints at the initial energy scale, quickly break one of these constraints once the RGE evolution is considered. This is also reflected in the plots in the lower row of figure 1, in which we show the points in the planes λ 1 -λ 3 on the left and (λ 4 − λ 5 )-(λ 1 + λ 2 + λ 3 ) on the right. In the left plot one can see that verifying the theoretical constraints exclusively at the initial scale µ = v gives rise to parameter points with values of λ 1 4 and −3 λ 3 8, whereas demanding that the constraints are respected within a range of energy of at least v ≤ µ ≤ 1 TeV, the allowed ranges shrink to λ 1 3 and −2 λ 3 5. 8 A similar observation can be made in the right plot, in which the allowed values change from λ 1 + λ 2 + λ 3 9 and −6 λ 4 − λ 5 10 to λ 1 + λ 2 + λ 3 6 and −4 λ 4 − λ 5 7 once the RGE running and the additional constraint µ v > 1 TeV is taken into account. Note that the limits on the values of the couplings that we found are much below the naive perturbativity criterion |λ i | < 4π, which is often applied in the analysis of extended Higgs sectors in order to exclude non-perturbative parameter regions.
Consequently, we conclude that regarding the collider phenomenology the main impact of our choice to demand the theoretical constraints to be respected at least until µ = 1 TeV is that the masses of the heavy states are closely related to the overall mass scale M , which, however, does not significantly constrain the values of λ 6,7,8 , since they do not depend directly on M (see appendix A). Thus, the only exception to the constraints on the mass splittings arises when there is a Higgs boson h 2 or h 3 with almost 100% singlet component present, in which case its mass would be dominantly related to the value of v S instead of M , and the mass could differ substantially from m A , m H ± and m H , as will also be further discussed below. Thus, our approach of including the theoretical constraints drives the model predictions towards the decoupling limit of the S2HDM, where the masses of the heavy states m A , m H ± and m H are approximately determined by the scale M of the soft-breaking of the discrete Z 2 . Considering the theoretical constraints described above has in some aspects the same effect as applying the constraints from the EWPO, which are also sensitive to large mass splittings between the scalar states [106]. This fact on its own is not very surprising since also the EWPO observables arise from the radiative corrections. More interesting, however, is that while it is sufficient to have either m H ∼ m H ± or m A ∼ m H ± in order to be in agreement with the constraints from EWPO (at one-loop level), the inclusion of the RGE running and the requirement µ v > 1 TeV gives rise to the fact that both conditions should be approximately fulfilled, i.e. m H ∼ m A ∼ m H ± .
The low values of µ v that we found for values of |λ i | 1 are relevant also for cosmological aspects of the S2HDM, where we stress again that one of the main motivations of the model is the possibility of accommodating a first-order EW phase transition. In order to achieve such a transition, it is required (just as in the 2HDM) to consider parameter space 8 λ1 has to be positive according to the BfB conditions on the tree-level scalar potential. regions where large loop corrections to the scalar potential are present, since at tree level the scalar potential does not allow for an EW phase transition of first order. The required loop corrections have their origin in values of one or more |λ 1,2,3,4,5 | > 1 [33]. As a result, our analysis indicates that for a perturbative study of the parameter regions of the S2HDM relevant for possible first-order EW phase transitions, it is of crucial importance to take into account constraints in relation to the perturbative unitarity and the RGE running of the quartic couplings. 9 On the other hand, if one restricts an analysis of the S2HDM to regions of the parameter space in which the couplings λ i have absolute values substantially below one, then the model can be valid to energy scales much beyond the TeV scale. In this case, however, the S2HDM cannot accommodate a first-order EW phase transition and its related phenomenology, and also the heavy BSM states are largely decoupled from the EW scale (as discussed above).
To shed more light on the spectrum of the Higgs bosons, we show in figure 2 the mass m h 3 in dependence of m h 2 on the left and m A in dependence on m H ± on the right. In the left plot one can see that it is possible that h 2 is substantially lighter than h 3 when it has a large singlet component of Σ h 2 > 90%, as indicated by the colours of the points. On the other hand, when h 2 and h 3 are sizably mixed, the masses of both states have to be relatively close to M in order to comply with the theoretical constrains. The same observation can be made in the right plot regarding the masses of A and H ± . Note here that all the points with m A − m H ± 150 GeV are grey, indicating that they feature values JHEP10(2021)215 of µ v < 1 TeV. In this plot the colour coding indicates the values of m h 2 , and a correlation can be seen between the mass of h 2 and the masses of A and H ± . The heavier the latter states, the larger also tend to be the values of m h 2 . Since by definition m h 2 < m h 3 , one can conclude that in most parameter points all six states h 2,3 , A and H ± are relatively close in mass, with the only exception being a very singlet-like state h 2 with m h 2 M , as mentioned earlier already. Hence, our analysis shows a trend towards the decoupling limit of the S2HDM, in which at low energies the model could become practically indistinguishable from the SM. In this case, the only possibility to observe a BSM effect would arise from the DM phenomenology or a possible invisible branching ratio of h 125 if the decay h 125 → χχ is kinematically allowed (see the discussion below). The presence of an invisible branching ratio of h 125 could also allow for a distinction between the S2HDM and the 2HDM, whereas the S2HDM in the decoupling limit could be practically indistinguishable from the pNG model with one Higgs doublet, in which case only a discovery of one of the additional particles of the S2HDM at a collider could shed light on the model realized in nature.
Going beyond the theoretical limitations, the spectrum of the Higgs bosons is also severely constrained by direct searches at colliders, where due to the fact that we focus here on the mass ordering with h 1 = h 125 only the LHC results play a role in our discussion. 10 Without going into the details of each of the relevant search channels, we list here the searches that were selected by HiggsBounds and which led to exclusions of parameter points in the scenario under investigation: Complementary to the direct searches for the BSM particles, the Higgs sector of the S2HDM can also be probed indirectly via the properties of the Higgs boson h 1 = h 125 resembling the Higgs boson that was discovered at the LHC. In order to illustrate the impact of such constraints, we show in the left plot of figure 3 the allowed parameter points, which all fulfill the criterion χ 2 125 ≤ χ 2 SM,125 + 5.99 (see section 3.2 for details), with sin(α − β) on the horizontal and tan β on the vertical axis. In the case in which one of the heavier states h 2 or h 3 has a singlet component of almost 100%, the S2HDM features an alignment limit similar to the 2HDM. In this limit the couplings of h 1 reduce to the ones of a SM Higgs boson, and the limit is determined by the condition sin(α − β) = 0 (see also ref. [76]). Consequently, departures from this condition are associated with deviations of the predictions for the signal rates of h 125 with respect to the SM. As can be seen in the left plot of figure 3, our analysis indicates that in order to be in agreement with the measured signal rates, one has to fulfill roughly | sin(α − β)| 0.1. The largest departures from zero are found for the lower end of the tan β range, whereas for larger values of tan β the allowed range of | sin(α − β)| shrinks substantially. The colour coding of the points indicates the singlet component of the SM-like Higgs boson Σ h 1 . Notably, we find that the current uncertainties of the signal-rate measurements still allow for a singlet-component of more than 14%.
Precise measurements of the properties of h 125 , for instance correlated deviations from the SM prediction of the various different couplings coefficients of h 125 to the up-and downtype fermions C h 125 uū and C h 125 dd and the gauge bosons C h 125 V V , could help to distinguish the type II S2HDM from the usual 2HDM. Here the coupling coefficients C h 125 uū,dd,V V are defined to be the couplings normalized to the ones of a SM Higgs boson. A sizable singletcomponent of h 125 , as found in parts of our parameter points, gives rise to a suppression of C h 125 V V . In the usual 2HDM, a deviation from |C h 125 V V | = 1 is possible via departures from the alignment limit, and thus tightly constrained to values of C 2 h 125 V V 0.9 [106]. Since we find parameter points with Σ h 125 > 0.1, and since in the S2HDM one has C 2 h 125 V V ≤ 1−Σ h 125 , a possible future measurement indicating C 2 h 125 V V 0.9 at the (HL)-LHC would favor an S2HDM interpretation instead of the 2HDM. It is also interesting to compare the maximum values of Σ h 125 ∼ 14% with the corresponding values found in the pNG DM model with only one Higgs doublet. In ref. [23] it was shown that in this case the mixing of the SM-like Higgs boson with the singlet state is more constrained, and, except when the singlet scalar and the doublet scalar are degenerate in mass, only values of up to 10% were found to be in The mixing among the CP-even scalar fields in the S2HDM is identical to the one of the N2HDM, such that it is not surprising that we find similar effects on the allowed parameter ranges of α i in the S2HDM. However, a crucial difference between both models is the presence of the additional particle χ in the S2HDM. Since we are focusing here on the Higgs funnel region of the model, it is possible that m χ < 125 GeV/2, giving rise to an additional decay mode of h 125 into an invisible final state. To illustrate the impact of this additional decay on the allowed parameter regions, we show in the right plot of figure 3 the branching ratio for the invisible decay of h 125 in dependence of m χ for the parameter points with m χ < 125 GeV/2. Here we show also the parameter points that would be excluded by the observed upper limit on BR(h 1 → χχ) [101] (indicated by the horizontal dashed line). In this way we can demonstrate the interplay between the global constraints from the HiggsSignals analysis and the direct limit on the invisible branching ratio. One can see that only a very small fraction of the otherwise allowed parameter points, which in -22 -

JHEP10(2021)215
particular have passed the constraint χ 2 125 ≤ χ 2 SM,125 + 5.99, lie above the ATLAS limit on the invisible branching ratio. Nevertheless, for some points we find values of BR(h 1 → χχ) that are about 50% larger than the upper limit in the whole range of m χ in which allowed points were found.
The grey points in the right plot of figure 3 have to be discarded because they feature a too large thermal relic abundance of DM. For the allowed points the DM relic abundance is indicated by the colour coding of the points. One can see that we find a limit of m χ ∼ 53.8 GeV below which no allowed points were found. 11 This limit arises from a combination of the upper limit on BR(h 1 → χχ), on the one hand, and the constraint Ωh 2 ≤ (Ωh 2 ) Planck , on the other hand. Parameter points with m χ 53.8 GeV feature either a χ that is weakly coupled to h 125 , in which case BR inv can be in agreement with the ATLAS limit but the DM relic abundance is too large because the annihiliation process with h 125 in the s-channel is not efficient (see also discussion below), or χ is coupled more strongly to h 125 , in which case the DM relic abundance can be below the upper limit but the invisible branching ratio of h 125 is unacceptably large. Note here that in the plot almost all parameter points with m χ 53.8 GeV belong to the first option, predicting too large values of Ωh 2 , while BR inv is below the experimental upper limit. On the other hand, there are only three parameter points with m χ 53.8 GeV belonging to the second option, featuring too large values of BR inv but with Ωh 2 below the Planck limit. The reason for this lies in our procedure to generate the parameter points using the genetic algorithm. Parameter points with values of BR inv = O(0.1) feature overall larger values of χ 2 125 and constitute therefore only a very small part of the sample of parameter points, because the genetic algorithm tries to find parameter points that minimize χ 2 125 (see the definition of the loss function defined in eq. (4.2)).
The above discussed findings already indicate the strong interplay between the Higgs phenomenology and the DM sector of the S2HDM, in particular in the scenario discussed here that fundamentally relies on the Higgs funnel to predict a DM relic abundance in agreement with experiments. To shed more light on this interplay, we show in figure 4 the relic abundance as predicted according to the freeze-out mechanism in dependence of the DM mass m χ . One can see the strong suppression of Ωh 2 for most parameter points at m χ ∼ 125/2 GeV, where the DM annihilation cross sections with h 1 in the s-channel are resonantly enhanced. At this precise resonance region, there are nevertheless also a few parameter points featuring values of Ωh 2 within an order of magnitude below the experimentally measured value (Ωh 2 ) Planck = 0.119 (indicated by the grey dashed line in figure 4). For these parameter points the resonant enhancement of the annihilation cross sections is counteracted by strongly suppressed couplings of χ to h 125 . 12 For values below this range the predicted amount of DM density is always too large (grey points). As already mentioned before, the reason for this lies in the constraints on the properties of h 125 . In order to predict an allowed value for Ωh 2 when m χ 53 GeV it is required that the coupling of χ to h 125 is large. However, this inevitably results in values of the invisible branching ratio for the decay h 125 → χχ above the experimental upper limit. As a result, the lower limit on m χ found here can be regarded as a robust bound under the assumption that h 125 corresponds to the lightest scalar h 1 . One can compare also to the right plot of figure 4, where the colour coding indicates the values of the singlet component of the SM-like Higgs boson h 1 . A clear distinction is visible between the points below and above the resonance at m χ = 125/2 GeV. Points with m χ below the resonance have substantially smaller values of Σ h 1 , whereas points with m χ above the resonance allow for values of Σ h 1 0.1. Moreover, only points for which m χ is relatively close to the kinematic threshold of the decay h 1 → χχ, i.e. m χ ∼ m h 125 /2, feature sizable values of Σ h 1 when m χ < m h 125 /2. The reason for this is that the couplings λ 7 and λ 8 that couple the singlet field to the doublet fields (see eq. (2.2)) appear in the partial decay width for the invisible decay as shown in eq. (3.1). In addition, these couplings are responsible for the possible singlet admixture of the state h 1 . Accordingly, parameter points with sizable values of Σ h 1 have sizable values of λ 7 and λ 8 , which in turn can give rise to too large values of BR(h 1 → χχ) whenever this decay is kinematically allowed. In section 4.2 we will address the question whether the bound m χ 53 GeV can be substantially modified in a scenario featuring a scalar h 1 with a mass smaller than 125 GeV, and the second lightest scalar h 2 plays the role of the discovered Higgs boson. In this case χ has two possibilities to annihilate resonantly, either with h 1 or with h 2 in the s-channel, and the predictions for the relic abundance can be substantially modified.

JHEP10(2021)215
For values of m χ > 125/2 GeV one can see that the prediction for Ωh 2 rises quickly with increasing value of m χ , because the resonant enhancement of the annihilation cross section is lost. As a result, most parameter points predict a too large DM relic abundance. Taking into account the values of m h 2 (indicated by the colour coding of the points in the left plot of figure 4), one can see that most parameter points with m χ 65 GeV that are in agreement with the upper limit on Ωh 2 feature a relatively light scalar h 2 with masses at the lower end of the scan range of m h 2 . As before, the reason for this is that when h 2 is not much heavier than twice the value of m χ , the second s-channel contribution to the annihilation cross section becomes relevant. This gives rise to a suppression of Ωh 2 such that the prediction can be below the experimental limit even when m χ is several GeV larger than 125/2 GeV. Again, this hints to the fact that also in the mass range m χ > 125/2 GeV the prediction for Ωh 2 could be substantially modified using the inverted mass hierarchy in which h 125 is not the lightest scalar, and we investigate this possibility assuming a Higgs boson h 1 at 96 GeV in section 4.2.
In both plots in figure 4 the grey points are characterized by either being excluded due to Ωh 2 > (Ωh 2 ) Planck , as already mentioned before, or they are excluded due to the constraints from DM indirect-detection experiments. In most parts of the analyzed parameter space, the more constraining experimental limit results to be the upper limit on the predicted relic abundance, as indicated by the fact that most of the grey points lie above the horizontal dashed line indicating the Planck measurement. However, there is a small region with 62.5 GeV m χ 67 GeV in which we find grey points below the Planck limit. Consequently, in this mass range of χ the indirect-detection limits from the observation of dSph by the Fermi satellite are more constraining. Note that this is a region in which it appears to be relatively easy to accommodate a value of Ωh 2 ∼ (Ωh 2 ) Planck without being in tension with constraints on h 125 , since it is just above the resonance of the annihilation cross section, and the decay h 125 → χχ is kinematically forbidden. The fact that these parameter points can be probed via indirect-detection experiments is therefore crucial. We remind the reader that the constraints derived from the Fermi measurements are subject to uncertainties, as was also discussed in section 3.2, such that the respective limits might change slightly in the future and are currently possibly not as robust as the Planck limit on the relic abundance. Nevertheless, our results indicate that when the DM candidate of the S2HDM in this mass range is responsible for a large fraction of the measured relic abundance, the observation of dSph and the resulting constraints (or signals, more optimistically speaking) will be of great importance for studies in the context of pNG DM.
In this context it is interesting to note that recent indirect-detection experiments found anomalies in the cosmic ray spectra coming from the center of the galaxy. The first so-called center-of-galaxy excess was found by the Fermi satellite, which measured an intensity of gamma-rays coming from the center of the galaxy significantly above the predictions of the standard model of cosmic rays generation and propagation with a peak in the spectrum around a few GeV [48,49]. Another anomalous cosmic-ray spectra was measured by the Alpha Magnetic Spectrometer (AMS) [62], mounted on the international space station, which reported an excess over the expected flux of cosmic ray antiprotons (see section 1 for details). 13 While it is still under debate whether the excesses arise from unresolved astrophysical sources [59][60][61] or the treatment of systematic uncertainties [136,137], or whether their origin could be the annihilation of DM, we will in the following assume that the latter is the case. 14 In ref. [66] it was shown that the excesses are compatible with a DM interpretation, where the DM candidate annihilates into pairs of b quarks. For the γ excess the allowed range of the mass of the DM candidate at the 2σ confidence level was found to be 37 GeV ≤ m DM ≤ 67 GeV. For thep excess the allowed range was found to be 46 GeV ≤ m DM ≤ 94 GeV, which partially agrees with the mass range preferred by the γ excess. Consequently, it is an interesting question whether the S2HDM can explain both the excesses simultaneously, while being in agreement with all theoretical and experimental constraints. 15 In order to answer this question we show in figure 5 on the vertical axes σ bb v rel , being the predicted velocity-averaged annihilation cross sections of χ into pairs of b quarks, in dependence of the DM mass m χ for the parameter points of our scan. The 2σ confidence level regions of these two parameters required to explain the γ and thep excesses are indicated with the blue and orange dashed lines, respectively [25,66]. We show the parameter points in the two plots of figure 5 under two different assumptions. In the left plot we assume that the usual thermal freeze-out scenario can be applied, such that we have to take into account the 13 The updated result of the AMS collaboration could neither definitively rule out nor confirm the DM interpretation of the antiproton excess [135].
14 See also refs. [138][139][140][141][142] for recent discussions of possible explanations of the center-of-galaxy excesses. 15 See ref. [25] for an investigation of the excesses in a singlet-extension of the SM featuring pNG DM.

JHEP10(2021)215
predicted values of the relic abundace for each parameter point. Hence, the values of σ bb v rel on the vertical axis are scaled by the factor ξ 2 as defined in eq. (3.2). On the other hand, in the right plot we show the parameter points under the assumption that the relic abundance of DM is always accounted for by χ, independently of the prediction from the thermal freeze-out. As a result, they demand a non-standard cosmological history giving rise to the experimentally measured relic abundance, which we will however not specify any further. In the plots the grey points correspond to parameter points that are excluded by a too large predicted relic abundance (left) or by constraints from dSph observations (left and right). In the right plot the dSph constraints are consequently applied also assuming Ωh 2 = (Ωh 2 ) Planck . Assuming the usual thermal freeze-out scenario (left plot), one can see that the resonant structure of the distribution of the annihilation cross sections gives rise to two distinct regions of m χ in which points inside the blue and the orange curves can be found. The first region at lower DM masses of m χ ≈ 50 GeV contains parameter points that predict values of ξ 2 σ bb v rel as required for an explanation of the excesses, and where the values of m χ lie roughly in the center of the values preferred by the γ excesses and at the lower end of the range preferred by thep excess. However, these points are excluded because the predicted values of Ωh 2 are about an order of magnitude larger than the experimentally measured value, as can also be seen in the left plot of figure 4. Accordingly, the parameter points in this region of m χ are excluded and the cosmic-ray excesses cannot be realized there. The second region of DM masses in which points within both the blue and the orange curves are found is given by 63 GeV m χ 67 GeV. However, as before, the corresponding points are shown in grey and are consequently excluded. Interestingly, here the responsible experimental constraint do not arise from the Planck measurement of the relic abundance, but from the Fermi-LAT observations of dSph, as was already discussed before. In fact, the predictions for Ωh 2 in this range of m χ are close or effectively identical to the Planck measurement. Hence, the points in this second region of DM masses possibly predict the correct DM relic abundance and could give rise to both the cosmic γ-and thep-excesses, but they are in tension with the null-results from the observations of dSph. Here we remind the reader, as was discussed already in section 3.2, that the Fermi-LAT dSph constraints are subject to uncertainties in regards to the astrophysical modelling of the spectral curves, and as a result might be slightly weaker as compared to applied here. Nevertheless, with future improvements of the dSph observations, for instance, due to the inclusion of more dSph and the increasing time periods of data taking, a firm exclusion (or confirmation if a DM signal will actually be found) of the parameter space region of the second DM mass range discussed here should be possible [143].
Under the assumption of a non-standard cosmological history that somehow gives rise to a relic abundance of χ in agreement with the Planck measurement (right plot), one can see that this time only one DM mass region with parameter points suitable for a realization of the excesses is present. Naturally, this region lies where the resonant enhancement of σ bb v rel is present, i.e. at 61 GeV m χ 67 GeV, which consequently partially coincides with the second region of DM mass found in the left plot of figure 5. As before the points that lie within both the blue and the orange curves are in tension with the dSph observations from the Fermi satellite.

JHEP10(2021)215
We end the discussion of the DM properties in this scan by noting that many of the above mentioned findings crucially depend on the assumed mass ordering of the CP-even Higgs bosons. In particular, the presence of a Higgs boson below 125 GeV can potentially impact the predictions for the relic abundance, as discussed in relation to figure 4. Moreover, the question whether the cosmic-ray excesses can be accommodated more easily when a second s-channel resonance for the annihilation cross section is available can be addressed. In section 4.2 we will investigate these questions following the approach of ref. [25], in which the presence of a Higgs boson at around 96 GeV was assumed in order to simultaneously explain also two collider excesses found at LEP in the bb final state and at the LHC in the diphoton final state.

pNG DM and a Higgs boson at 96 GeV
In ref. [25] it was used that the hypothetical particle state h 96 at 96 GeV can be coupled to new relatively light charged states that can give rise to additional contributions to the loop induced coupling of h 96 to photons in order to account for the diphoton excess found by CMS. In ref. [69] it was shown that in the N2HDM the presence of the additional doublet Higgs field and the real singlet field are sufficient to accurately describe the collider excesses. Here, the diphoton rate was enhanced not via an enhancement of the coupling coefficient |C h 96 γγ |, where the coupling coefficients C h 96 ... are defined as the coupplings normalized to the one of a SM Higgs boson of the same mass. Instead, the branching ratio for the diphoton decay of h 96 was enhanced via a suppression of the couplings of h 96 to b quarks, which then also gives rise to a suppression of the total width of h 96 . 16 The required suppression of the coupling coefficient |C h 96 bb | (without suppressing |C h 96 tt | in order to maintain sizable couplings to photons via the t-quark loop) can also be realized in the S2HDM due to the possible mixing patterns in the CP-even sector and the presence of the three mixing angles α 1,2,3 in total analogy to the N2HDM. In this regard, the only difference in the S2HDM compared to the N2HDM is the possible presence of the additional decay modes h 96 /h 125 → χχ, potentially giving rise to a small suppression of the decay modes h 96 → γγ relevant for the CMS excess and h 96 → bb relevant for the LEP excess, or to stronger constrains on the properties of h 125 . In the following we will discuss a scan to illustrate the impact of the presence of h 96 on the phenomenology of the DM candidate χ, and whether the collider excesses can be realized in combination with the cosmic-ray excesses.
Before going into the description of the parameter scan that we performed, we briefly introduce the relevant details of the collider excesses. At LEP searches for Higgs bosons were performed utilizing the bb final state [144], which can be exploited at a lepton collider in contrast to the LHC due to the much smaller SM background. Theoretically, the Higgs boson that is searched for is assumed to be produced via the Higgstrahlung pocess and subsequently decays into a pair of b quarks. A local excess of about 2σ confidence level was observed at a mass of roughly 96 GeV, where the mass resolution is rather poor due to the hadronic final state. In ref. [70] it was shown that the excess is consistent with a signal 16 The presence of a second doublet field gives rise to the presence of the states H ± , such that also in the S2HDM (compared to the SM) new charged states are present. However, the loop contributions of H ± to |C h 96 γγ | are not relevant for the explanation of the CMS excess, such that one can have m H ± 96 GeV. Low-mass Higgs-boson searches have also been performed at the LHC in various final states. CMS searched for light Higgs bosons in the diphoton final state utilizing the 8 TeV and parts of the 13 TeV datasets [145]. A local excess of roughly 3σ confidence level was observed at a mass of 96 GeV, hence in agreement with the mass range compatible with the LEP excess.
In this case the excess is consistent with a signal interpretation corresponding to a signal strength of In our scan, in which h 1 will play the role of the state h 96 , we compare the theoretical predictions for the signal strengths to the experimental values given above. The predictions were calculated by Hence, in both cases the cross section ratios that enter the definitions of the signal strengths are expressed to a very good approximation in terms of the effective coupling coefficients in order to quantify the goodness of the fits to the excesses. In this definition we assumed that there is no correlation between both measurements. Technically, the details of the scan that we discuss here are very similar to the ones of the scan discussed in section 4.1. The scan ranges were set as given in eq. (4.1), except for the masses of the scalars, which were chosen to be such that m h 3 = m H is further constrained by the condition ∆M max < 200 GeV, as defined in eq. (4.1) and substantially heavier than h 1 and h 2 due to the lower limit on m H ± . We again followed the two-step procedure. In the first step, we used the genetic algorithm to obtain parameter points in agreement with the theoretical constraints and the experimental constraints from the Higgs phenomenology. To the loss function defined in eq. (4.2) we added a term 10χ 2 96 in order to obtain parameter points that potentially feature both a good fit to the signal rates of the SM-like Higgs boson h 2 = h 125 and to the signal rates  where the value of χ 2 SM,96 is obtained from eq. (4.6) assuming zero values for both µ LEP and µ CMS as predicted by the SM, in which no particle is present at a mass of 96 GeV. As a result, in comparison to the analysis discussed in section 4.1 in which the requirement χ 2 125 ≤ χ 2 SM,125 was used, the requirement shown in eq. Here it should be noted that even in the most extreme case with χ 2 96 = 0 the allowed maximum value of χ 2 125 still does not indicate severe modifications of the signal rates of h 125 , taking into account that the HiggsSignals fit result applies a total amount of 107 observables, such that the reduced χ 2 value remains substantially smaller than one even in this case. The second step is totally analogue to the scan discussed in section 4.1. All parameter points that pass the constraint shown in eq. (4.8) were confronted with the theoretical constraints including now the RGE evolution of the parameters. As before, we required the scalar potential to be well behaved up to energy scale of at least µ v = 1 TeV, such that in particular the values of the quartic couplings λ i allow for a perturbative treatment at the range of energy at which there are also particle masses in our scan. Finally, the remaining experimental constraints regarding the DM phenomenology were applied.
We show the resulting parameter points in figure 6, where we display the signal rate µ LEP in dependence of µ CMS . We indicate with the colour coding of the points the value of the DM mass m χ (left) and the DM relic abundance Ωh 2 as predicted by the usual thermal JHEP10(2021)215 freeze-out scenario (right). Also shown as grey points are parameter points that are excluded by a too large prediction of the relic abundance or by limits coming from observations of dSph. The ellipse in both plots indicates the region in agreement with the collider excesses at the 1σ confidence level, i.e. χ 2 96 = 2.3. One can see that we find parameter points within the ellipses. Consequently, both excesses can be explained simultaneously while taking into account the constraints described in section 3. In the left plot, we observe that parameter points with sizable values of µ LEP and µ CMS feature DM mass values close to or larger than m h 125 /2. On the other hand, parameter points with m χ < m h 125 /2 only predict substantially smaller signal strengths, and the collider excesses cannot be accounted for. The reason for this is, as was also discussed in section 4.1, that in this case the decay h 125 → χχ is kinematically open. As a result, the possible mixing of the singlet field h 1 = h 96 with the SM-like Higgs boson h 2 = h 125 is much more constrained. However, a sizable mixing of h 96 and h 125 is necessary to obtain values of µ LEP and µ CMS of the order of the experimentally measured values. We therefore can conclude that a realization of the collider excesses demands DM masses of m χ > m h 125 /2. In the right plot of figure 6 we find that several of the parameter points that are able to explain both excesses also predict sizable values for the relic abundance, with some parameter points saturating the value measured by the Planck collaboration. Accordingly, we come to the conclusion that the S2HDM can accommodate the collider excesses at 96 GeV while at the same time accommodating a large fraction or all of the measured DM relic abundance.
In figure 7 the predicted relic abundance is shown in dependence of the DM mass. The values of the signal rates measured by LEP (left) and CMS (right) are also indicated by the colour coding of the points. We note a new prominent feature in the distribution of the parameter points with respect to figure 4. Due to the opening of a new resonant s-channel mediated by the h 96 , parameter points featuring DM masses smaller than about 53 GeV can now be in agreement with the upper limit imposed by the observed DM relic abundance. Moreover, the presence of h 96 also gives rise to the fact that a large fraction of parameter points with m χ > m h 125 /2 lie below the Planck limit, whereas we found in section 4.1 (compare to figure 4) that in this DM mass region most points predict Ωh 2 > (Ωh 2 ) Planck . Grey points that lay below the experimental upper limit are excluded by dSph observations. Here it is interesting to note that we find, in addition to the region around m χ ∼ 63 GeV already present in figure 4, a second region at 48 GeV m χ 58 GeV in which the dSph constraints discard points that would be in agreement with the Planck measurement of the DM relic abundance. In the left plot of figure 7 we find that for the points at the right side of the resonance the predicted values of µ LEP can be close to the measured central value µ exp LEP = 0.117 independently of the precise value of m χ . On the contrary, as can be seen in the right plot of figure 7, values of µ CMS ∼ µ exp CMS = 0.6 that are in agreement with the constraints are mostly found in the interval 62 GeV m χ 65 GeV. For larger values of m χ one can still find parameter points that fit the CMS excess at the level of 1σ. However, they often predict too large values of Ωh 2 > (Ωh 2 ) Planck and are therefore shown mostly as grey points. The reason for this is that, as discussed before, fitting the diphoton excess requires a suppression of the couplings of h 96 to b quarks. However, this then yields also a suppression of the annihilation cross section via the process χχ → h 96 → bb.  In order to discuss the gamma-ray and the antiproton excesses, we show in figure 8 today's velocity-averaged annihilation cross section of χ into pairs of b quarks taking into account the number density as predicted by thermal freeze-out (left) and assuming Ωh 2 = 0.12 (right), as explained in section 4.1. In comparison to figure 5, here we observe -32 -

JHEP10(2021)215
that there are more regions of m χ in which points are found inside the preferred region to explain both cosmic-ray excesses simultaneously. These points remain in tension with present limits imposed by the observation of dSph. We remind the reader about the uncertainties in determining those limits (see section 3.2 for more details). Regarding the agreement with the signal rate µ CMS , only the parameter points situated towards the right end of the blue curve could simultaneously explain the two cosmic ray and the CMS excesses. These points are again in tension with indirect-detection limits from dSph observations. Regardless of whether the collider excesses are accommodated or not, we see that the presence of h 96 gives rise to more points at the lower end of m χ that lie within the blue and the orange curves. Thus, the new light scalar state gives rise to new interesting regions of parameter space with m χ < 60 GeV in the context of the cosmic-ray anomalies. However, as was already mentioned, the collider excesses, which were the main motivation to investigate a scenario with m h 1 = 96 GeV in the first place, cannot be realized here.

Conclusions
In this paper we analyzed a singlet-extended 2HDM, called S2HDM, which is a model with a rich Higgs phenomenology and that incorporates a pseudo-Nambu-Goldstone boson dark matter candidate. We focused on the parameter space of the S2HDM featuring DM masses in the Higgs funnel region, i.e. 40 GeV ≤ m χ ≤ 80 GeV. One of the main purposes of this analysis was to illustrate the combined impact of various theoretical and experimental constraints on the model parameters, where in particular the strong interplay between the Higgs-sector phenomenology and the DM sector of the S2HDM was demonstrated. We required the scalar potential to be well-behaved up to energy scales of 1 TeV, i.e. to be bounded-from-below, to feature a stable electroweak vacuum and to fulfill conditions derived from perturbative unitarity. We also ensured that the parameter points were in agreement with measurements of electroweak precision observables, flavour physics, properties of the discovered Higgs boson at 125 GeV, searches of additional scalar states and the DM observables. The model exploration of the multi-dimensional parameter space of the S2HDM was performed with the help of a genetic algorithm, by which, compared to random scans of the parameters, a significant improvement of the computing time required to find viable parameter points was achieved.
In our numerical analysis, we focused on two benchmark scenarios. Firstly, we performed a broad parameter scan assuming that the SM-like Higgs boson h 125 was the lightest of the three CP-even Higgs bosons, such that the predictions for the DM relic abundance assuming the usual thermal freeze-out mechanism are mainly determined by the resonant s-channel annihilation mediated by h 125 . Secondly, we studied a scenario featuring a singlet-like CP-even state h 96 at 96 GeV, where the presence of h 96 also gives rise to a second s-channel contribution to the thermal freeze-out cross section and today's annihilation cross section relevant for DM indirect-detection experiments.
In the first scenario in which h 125 is assumed to be the lightest CP-even scalar, DM masses 62.5 GeV m χ 67 GeV were found to be able to explain the γ-ray and antiproton cosmic rays excesses, while simultaneously also predicting values of the DM relic abundance -33 -

JHEP10(2021)215
in agreement with the observations by the Planck collaboration. However, these parameter points are in tension with indirect-detection limits derived from observations of dwarf spheroidal galaxies, where it should be taken into account that these constraints are still subject to uncertainties in regards to the astrophysical modelling of the spectral curves. Regarding the Higgs phenomenology, we found that demanding that the theory can be treated perturbatively up to energy scales of at least 1 TeV has a strong impact on the Higgs spectra that can be realized. Namely, the mass splittings among the heavy scalar states H ± , A and h 2,3 were found to be smaller than roughly 100 GeV, driving the model towards the decoupling limit, with the only exception of a very singlet-like CP-even state h i , which can be substantially lighter (or heavier) than the other BSM state without giving rise to issues with unperturbative effects at energy scales below 1 TeV.
In the second scenario, we studied whether the S2HDM could offer an explanation for the collider excesses observed at about 96 GeV at LEP and CMS in the bb and the diphoton final state, respectively. Here we found that a singlet-like CP-even Higgs boson at 96 GeV can reproduce both collider excesses under the constraint that m χ > m h 125 /2 in order to allow for a sizable mixing between h 96 and h 125 . Furthermore, it is possible to accommodate at the same time a large fraction or all of the measured DM relic abundance. Finally, we found that the simultaneous explanation of the cosmic-ray excesses and the collider excess at 96 GeV is in principle possible, but, as in the first scenario, the parameter points are also in tension with limits arising from observations of dwarf spheroidal galaxies.
To summarize, we demonstrated that the S2HDM is an attractive model that can accommodate a rich phenomenology and an interesting interplay between the DM sector and the Higgs sector. We also showed that it is crucial to take into account the various theoretical and experimental constraints on the model parameters. For future studies, we make our implementation of the model predictions and the application of the constraints available to the public in the form of a python package called s2hdmTools, which is briefly described in appendix C.

C The python package s2hdmTools
In order to make the analysis of the S2HDM as performed here publicly available, we developed the python package s2hdmTools. The code can be downloaded at https://gitlab. com/thomas.biekoetter/s2hdmtools. The installation requires a python3 environment and compilers for Fortran, C and C++. In addition, a python2 installation is required for the external code MadDM. During the installation process, the following external codes are downloaded and installed: -

JHEP10(2021)215
The user interface of s2hdmTools is defined in the class ParamPoint and briefly summarized in the README. The most important features can be accessed via the following functions defined in ParamPoint: More instructions regarding the installation and the usage of the package can be found in the README of the repository. In addition, the application folder contains several of the scripts that were used in order to produce the results discussed here.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.