Global fit of pseudo-Nambu-Goldstone Dark Matter

We perform a global fit within the pseudo-Nambu-Goldstone Dark Matter (DM) model emerging from an additional complex scalar singlet with a softly broken global U (1) symmetry. Leading to a momentum-suppressed DM-nucleon cross section at tree level, the model provides a natural explanation for the null results from direct detection experiments. Our global fit combines constraints from perturbative unitarity, DM relic abundance, Higgs invisible decay, electroweak precision observables and latest Higgs searches at colliders. The results are presented in both frequentist and Bayesian statisical frameworks. Furthermore, post-processing our samples, we include the likelihood from gamma-ray observations of Fermi -LAT dwarf spheroidal galaxies and compute the one-loop DM-nucleon cross section. We find two favoured regions characterised by their dominant annihilation channel: the Higgs funnel and annihilation into Higgs pairs. Both are compatible with current Fermi -LAT observations, and furthermore, can fit the slight excess observed in four dwarfs in a mass range between about 30–300 GeV. While the former region is hard to probe experimentally, the latter can partly be tested by current observations of cosmic-ray antiprotons as well as future gamma-ray observations.


Introduction
The true particle nature of Dark Matter (DM) continues to remain a mystery despite a plethora of astrophysical/cosmological evidence to support its existence [1]. Although the well-known Weakly Interacting Massive Particle (WIMP) offers a viable solution, most common models are strongly constrained by direct detection experiments [2][3][4]. This has forced us to either seek alternate particle DM candidates (e.g., axions [5][6][7][8], sterile neutrinos [9,10]) or explore new ways of saving the canonical 'WIMP paradigm'.
A natural way of achieving the latter is to suppress the DM-nucleon interaction at tree level. For instance, in certain particle DM models, some parameter combinations can lead to blind spots in direct detection experiments or even a suppression of the DMnucleon coupling [11][12][13][14]. Alternatively, the DM-nucleon couplings could vanish due to -1 -

JHEP04(2020)015
into HH and hH (where h and H are the 125 GeV Higgs and new scalar, respectively) as well as those proceeding via off-shell vector bosons. In addition, we compute the oneloop pNG DM-nucleon cross section for our samples and compare the resulting values against the current limits from XENON1T (2018), and projected future limits from LZ and DARWIN. Our FeynRules [44], UFO [45], CalcHEP [46] and FeynArts [47] model files are publicly available at the FeynRules database. 2 The remainder of the paper is organised as follows. In section 2, we introduce the pNG DM model. In section 3, we describe the various observables and likelihoods used in our global fit and in the post-processing. Our numerical scan details, global fit results, including Fermi -LAT and direct detection constraints, are presented in section 4. We conclude in section 5. Appendices A and B summarise analytic expressions used in this paper.

Pseudo-Nambu-Goldstone Dark Matter
We extend the SM Lagrangian by adding a new complex scalar field S that couples to the SM particles via a Higgs portal term, Φ † Φ (Φ is the SM Higgs doublet). The model Lagrangian is given by [26] L = L SM + L S + L soft , (2.1) where L SM is the SM Lagrangian, S → e iα S, (2.4) where α is a real, space-time independent parameter. However, the µ 2 S term in eq. (2.3) softly breaks this symmetry. Thus, the model contains a massive Goldstone boson, i.e., a pNG boson. After this symmetry breaking, we are left with a residual Z 2 symmetry, S → −S, of the dark U(1) group, which forbids a linear term in S in the above Lagrangians.
The parameter µ 2 S can be made real and positive by the phase redefinition of S. Thus, eq. (2.1) is invariant under a dark CP symmetry: This symmetry is unbroken by the S vacuum expectation value (VEV) as for positive µ 2 S , the VEV is real. Thus, the total symmetry of the model Lagrangian is Z 2 ⊗ CP .
With an extra scalar, the scalar potential becomes

JHEP04(2020)015
where V S and V soft can be read directly from eqs. (2.2) and (2.3) respectively. Meanwhile, the SM part of the potential reads (2.7) After electroweak symmetry breaking, the model spectrum can be analysed by decomposing Φ and S in the unitary gauge as where v h is the SM Higgs VEV. Under the dark CP symmetry in eq. (2.5), χ → −χ. This guarantees the stability of χ and makes it a viable DM candidate; the physical χ mass is m 2 χ = µ 2 S . After imposing the stationary point conditions at (φ, s) = (0, 0), we get Given that the S VEV is non-zero in general, the λ ΦS term in eq. (2.2) leads to a mixing between the CP -even interaction eigenstates (φ, s). Thus, the squared mass matrix M 2 is non-diagonal, namely As M 2 is real and symmetric, it can be diagonalised by a unitary transformation: where O = cos θ sin θ − sin θ cos θ . (2.13) Here θ is a mixing angle that satisfies the following relation: (2.14) The eigenvalues of M 2 correspond to the masses of the CP -even mass eigenstates (h, H), Given the discovery of a SM-like Higgs at the LHC [48,49] The remaining parameters in eqs. (2.2) and (2.3) can be expressed as (2.20)

Observables and constraints
In this section, we describe the set of constraints included in our global fit and postprocessing of final samples.

Theoretical bounds
We require the model parameters to satisfy the following two theoretical bounds.

1.
Bounded tree-level potential : the tree-level potential in eq. (2.6) must be bounded from below. This translates into the following lower bounds: 2. Perturbative unitarity: we require the perturbative unitarity of scattering amplitudes [50]. Using the HH → HH scattering process, we impose the following upper bound on the S quartic coupling [51]: Although this bound can vary with the exact scattering process of interest, we choose this form to maintain comparability with previous studies in literature [26,28,31,41].
Parameter points that do not fulfill these requirements are discarded from our scan. This is formally achieved by assigning a very small likelihood to such points.

Thermal relic abundance
The pNG boson χ is the DM candidate. Similar to the extended scalar singlet model [52], χ can annihilate into f f (where f = quarks/leptons), W + W − , ZZ, hh, hH and HH final states via an s-channel h/H exchange. In addition, χ annihilation into hh, hH and HH final states is also possible via t-and u-channels via χ exchange.

JHEP04(2020)015
In our numerical scans, we require χ to make up all of the observed DM relic abundance. 3 This is achieved using a Gaussian likelihood function for the DM relic density that is centered at the Planck (2018) measured value [53]: We also include a 5% theoretical uncertainty and combine it in quadrature with the Planck measured uncertainty. This is done to reflect any uncertainties arising from the relic density calculation in micrOMEGAs v5.0.8 [54].
With two neutral scalar mediators, the χ annihilation cross section is resonantly enhanced when m χ ∼ m h, H /2. To obtain the correct DM abundance, the χ annihilation cross section must be sufficiently suppressed. This is achieved for small values of v h /v s (or large v s ); an expression for the DM-scalar coupling can be found in appendix A. Away from these resonances, large values of v h /v s (or small v s ) generally saturates the χ relic density to the observed value.

Higgs invisible decay width
When m χ m h, H /2, the two scalars {h, H} are kinematically allowed to decay into a pair of DM particles, i.e., h, H → χχ. This contributes to the following invisible decay widths [28]: Recently, both the ATLAS [55] and CMS [56] experiments released new upper limits on the Higgs invisible branching ratio BR(h → χχ) for a SM-like Higgs from a combination of Run 1 and 2 analyses. Here we adopt the conservative upper limit from the ATLAS experiment [55], namely where Γ tot h (m h ) is the total decay width of h into SM and non-SM final states. In the following, we apply the limit in eq. (3.6) only on the scalar h whose mass is fixed at 125 GeV. This experimental limit is derived from Higgs production in association with a weak gauge boson or through vector boson fusion (VBF), and assuming a SM-like Higgs except for the fact that it can decay to a pair of invisible particles, e.g., DM. The experimental signatures are large missing energy with either a weak boson or a pair of jets. Thus, the invariant mass of the invisible particles is not measured, and the second JHEP04(2020)015 scalar can contribute as well to these processes. However, this contribution is small when the mixing angle is small as the production of the second (first) scalar is suppressed by a factor sin 2 θ (cos 2 θ) compared to the SM production rate. The production rate also varies with the mass of the scalar but this effect is small for a light scalar (m H ∼ 100 GeV) in the most constraining channel, i.e., VBF, as the cuts on the invariant mass of the two jets (m jj > 1 TeV) are already requiring a very large partonic center-of-mass energy. As the second scalar mass increases, its production rate is further reduced. To sum up, our approximation is only valid for small mixing angles and for m H = m h . However, as we show in the results section, the mixing angle is allowed to be large and even maximal when the two scalars are degenerate, i.e., when m H ∼ m h = 125 GeV. In this case, both scalars contribute to the process and interfere quite strongly. The amplitude can be written as where Γ tot h (Γ tot H ) are the total decay width of scalars h (H). The amplitude for the production of the scalar is multiplied by a factor cos θ (sin θ) compared to the SM due to the modification of the couplings between the gauge bosons and h (H). The remaining factors, besides the propagators, are due to the couplings with the DM particle χ (see appendix A). The missing pre-factor in eq. (3.7) depends only on pure SM couplings and its exact expression depends on the process considered. The two terms in the brackets can cancel exactly for θ = π/4 if the two masses are identical, as in that case, the two widths are also identical. Thus, these points are unconstrained experimentally but many of them would be excluded by applying blindly the constraints on the invisible decay width, which is the dominant channel for low values of m χ and v s . However, the correct re-interpretation of the invisible width constraint goes beyond the scope of this paper and thus is left as a future work.
The upper limit in eq. (3.6) constrains the m χ m h /2 region where Γ inv (h → χχ) can be sizeable. In our numerical scans, we use a one-sided Gaussian likelihood function that is centered at the above measured value for BR(h → χχ). Similar to the relic density likelihood, we add a 5% theoretical uncertainty from our calculation of BR(h → χχ) and combine it in quadrature with the (expected) branching ratio uncertainty of 0.07 [55].

Electroweak precision observables
The extra scalar S contributes to the gauge boson self-energy diagrams. Its effect on the electroweak precision observables (EWPO) can be parametrised by the oblique parameters S, T and U [57,58]. As S is electrically neutral, only the W and Z boson self-energies are modified.
Given that only the real component of S acquires a non-zero VEV, the pNG DM χ (imaginary component of S) does not contribute to the W/Z self-energies. Thus, the oblique parameters in our model have the same functional dependency as in the extended scalar singlet model [52], namely,  [60] ∆S = 0.04 ± 0.11, ∆T = 0.09 ± 0.14, ∆U = −0.02 ± 0.11, (3.9) and the following correlation matrix: (3.10) In our numerical scans, we use the following EWPO likelihood function [61]: where ∆O i are the central values for the shifts in eq. (3.9), Σ 2 ij ≡ σ i ρ ij σ j is the covariance matrix, ρ ij is the correlation matrix in eq. (3.10) and σ i are the associated errors in eq. (3.9).

Higgs searches at colliders
In the narrow-width approximation, the signal strength µ h for a SM-like Higgs h [28] is (3.12) The inclusion of H in our model leads to a universal suppression of couplings between h and SM particles. Thus, the h production cross section is where σ SM pp→h (m h ) is the h production cross section in the SM. Similarly, the branching ratio of h into SM particles is , (3.14) where Γ SM h (m h ) is the total decay width of a SM-like Higgs h with mass m h into SM final states. The last two terms in the denominator correspond to the new decay modes of h, namely h → χχ and h → HH.
Thus, the signal strength µ h in eq. (3.12) becomes is the h signal strength in the SM. From the above expression, it is clear that µ h = µ SM h when the mixing angle θ = 0, or when h decay into non-SM final states is kinematically allowed.
To constrain the parameters of the pNG DM model in our scan, we take into account the following two contributions to the likelihood:

JHEP04(2020)015
• We consider Higgs searches performed at the LEP experiment by utilising the Hig-gsBounds v5.3.2beta [62,63] package. It allows us to compute a χ 2 LEP for the mostsensitive LEP analysis on the basis of the scalar masses, effective Higgs-SM couplings, total decay widths and branching ratios. The corresponding likelihood function reads: • We take into account constraints from the observed Higgs signal strengths and mass measurements performed for a SM-like Higgs at the LHC. This is achieved using the HiggsSignals v2.2.3beta [64] package. In practice, we compute three contributions to the χ 2 that are based on i ) combined run 1 results (χ 2 R1 ); ii ) results from 13 TeV LHC analyses (χ 2 13 TeV ); and iii ) results in the form of Simplified Template Cross Sections (STXS) (χ 2 STXS ). The final HiggsSignals likelihood that we use in our scans is Each of the individual chi-squares are computed using the peak-centered method with Gaussian probability density function and zero theoretical mass uncertainty for the two scalar masses. For more details, see ref. [64].
The LEP analyses that we use are sensitive to parameter points with m H 120 GeV only, whereas the observed Higgs signal strengths potentially constrain points with m H in the whole considered mass range. 4

Fermi -LAT gamma-ray observations
Gamma-ray observations of dSphs by Fermi -LAT provide robust limits on the DM annihilation flux. To constrain our model, we use the publicly available energy-binned likelihood profiles 5 from Fermi -LAT [34], as implemented in MadDM v3.0 [65]. We consider the set of all 45 observed dSphs and use the measured J-factors based on spectroscopic observations [34] as adopted from ref. [66]. When measurements are not available, we use the values predicted from the distance scaling relationship with a nominal uncertainty of 0.6 dex [34]. We profile over the J-factor of each dwarf galaxy according to its uncertainty and obtain a total likelihood function as described in ref. [67]. The corresponding gamma-ray energy spectra are computed using MadDM, while showering and hadronisation is achieved using Pythia v8.0 [68]. We include all annihilation channels that contribute at least 1% to the total DM annihilation rate today. In particular, besides the usual 2 → 2 processes in MadDM, we also include 2 → 3 annihilation processes, χχ → V V * , where V * is an off-shell weak gauge boson. In addition, we include 2 → 3, 4 Note that we do not consider constraints from direct searches for a second (heavy) Higgs performed at the LHC (also implemented in HiggsBounds v5.3.2beta). Besides the technical limitation that HiggsBounds does not provide a likelihood for these searches, we found that the respective 95% C.L. limits are weaker than those of other constraints applied in our analysis. In particular, the EWPO constraint excludes sizeable values of θ in the region mH 130 GeV. 5 https://www-glast.stanford.edu/pub_data/1203/ -9 -JHEP04(2020)015 2 → 4, 2 → 5 and 2 → 6 processes such as χχ → Hh, HH where H decays further into pairs of SM particles, including V V * . The decay of all on-shell SM particles is performed within Pythia.
In four (Reticulum II, Tucana III, Tucana IV and Indus II) of the 45 dSphs, slight excesses have been found with a local significance of roughly 2σ each [34,69,70]. Consequently, a combination of the likelihoods from all dwarfs favors a certain range of DM masses and annihilation cross sections that provide a flux compatible with the excess. However, as the DM origin of these excesses is not yet established, we show our results in section 4 with and without including the respective four dwarfs. The latter choice imposes an upper bound on the DM annihilation cross section only.

Direct detection at one-loop level
The elastic spin-independent (SI) pNG DM-nucleon cross section is momentum-suppressed at tree-level (see appendix A). It is given by [30] σ tree-level Higgs-nucleon coupling [71][72][73], m N = 939 MeV is the averaged nucleon mass and v χ is the DM velocity in the laboratory frame. In the vicinity of the Earth, v χ ∼ 10 −3 . Thus, the nuclear recoil rate is suppressed by a factor of v 4 χ ∼ 10 −13 . For a typical choice of model parameters, the tree-level cross section in eq. (3.18) is too small to be experimentally observed at current or future planned experiments [30].
The first non-vanishing contribution to the DM-nucleon cross section at zero velocity appears at one-loop level. It can be approximated by [26,30] (3. 19) In the limit of m 2 χ ≡ µ 2 S → 0, the DM particle χ becomes a true Goldstone boson of the spontaneously broken global U(1) symmetry. Thus, the direct detection amplitude should vanish. This behavior is observed in eq. (3.19).
In ref. [30], the authors point out that the approximate one-loop cross section in eq. (3.19) can under/overestimate the actual cross section by several orders of magnitude depending on the model parameters. Thus, we use the full one-loop cross section from ref. [30] in our study: -10 -

JHEP04(2020)015
where the one-loop function F is Here the C and D terms are Passarino-Veltman functions [74][75][76] which we compute using LoopTools v2.14 [76,77]. 6 The coefficients A i are defined as Note that F is proportional to m 2 χ , and the fact that both the C 2 and D 3 functions behave as constants in the limit of m χ → 0 [26,30], the direct detection amplitude indeed vanishes in this limit.
By means of our global fit, we are able to check the conclusions of ref.
[30] more generally. We post-process our final samples, compute the approximate and actual oneloop cross sections using eqs. (3.19) and (3.20) respectively, and present results in the {m χ , σ 1-loop χN }-plane. This allows us to confront the allowed parameter space of the model against the current limits from XENON1T [4], and projected sensitivities from LZ [32] and DARWIN [33]. These results are presented in section 4.
The authors of ref. [31] have also computed the full one-loop DM-nucleon scattering cross section. Their results are consistent with ref. [30] in most parts of the parameter space. However, large deviations appear at small DM masses. We have verified that for the parameter space that we consider, these deviations do not impact our conclusions. In regions where these differences can be sizeable, the overall cross section lies well below the experimental limits that we show in section 4.

Results
In this section, we discuss the statistical treatment of the various constraints included in our global fit and show the corresponding results.

Statistical analysis
To find the allowed regions in the model parameter space, we use MultiNest v3. 10 tool designed to compute the Bayesian evidence Z (defined below). As a by-product, it draws posterior samples from a distribution that may contain a high multiplicity of nodes and/or degeneracies. Furthermore, MultiNest is capable of sampling the profile likelihood ratio (defined below) for the purpose of frequentist analysis.
A key ingredient for both a frequentist and Bayesian analysis is the likelihood (or log-likelihood) function. In our numerical scans, the total log-likelihood function is where θ ≡ (m χ , v s , θ, m H ) are the free model parameters. Each of the individual likelihood functions are described in section 3.
The range and prior types for our free model parameters are summarised in table 1. For the mixing angle θ, we find that our results are symmetric under θ → −θ. In addition, the case θ = π is analogous to θ = 0, thus we only restrict θ ∈ [0, π/2]. To cover the region close to the two resonances, m χ m h,H /2, where the annihilation cross section is enhanced, we scan up to very large values for the second scalar VEV, i.e., v s ∈ [10, 10 6 ] GeV. The upper boundary corresponds to very small couplings λ S and λ ΦS , see eqs. (2.19) and (2.20) respectively.

Profile likelihoods
In a frequentist analysis, the statistical precision of a parameter estimate is represented by a confidence interval that encapsulates the frequentist 'coverage probability'. Such an interval is dependent on the data x, and thus changes upon each re-iteration of the experiment. As proper frequentist coverage is usually not possible for complicated likelihoods and parameter spaces, approximate methods are often used [79]. One such method is the well-known profile construction [80], which depends on the profile likelihood ratio (PLR): Hereν(θ i , θ j ) are the parameter values {θ k | k = i, j} that maximise L total (θ) for a fixed (θ i , θ j ), whereasθ is the maximum likelihood estimate for θ, i.e., a 'best-fit' point that maximises L total (θ) [81,82]. To construct confidence intervals, we maximise Λ in the relevant parameter planes of interest while profiling over the other parameters and construct iso-likelihood contours at fixed confidence level (CL), e.g., 68 Observables Table 2. A summary of the best-fit (BF) points, key DM observables (the DM relic abundance, the dominant annihilation channel during freeze-out (FO) and today, the DM annihilation cross section today and the one-loop DM-nucleon cross section) and total log-likelihood ln L BF total (θ) from our global fit (column 1 ), and after post-processing our samples with Fermi -LAT likelihood with 41 (column 2 ) and 45 (column 3 ) dwarf spheroidal (dSph) galaxies.
In figure 1, we show our PLR plots in six 2D planes spanned by all combinations of four model parameters. These are generated using pippi v2.0 [83]. In each plane, model parameters that are not shown are profiled over. The 1σ (2σ) CL contours are marked by solid (dashed) lines. The best-fit point is shown as a red star; it is also summarised in column 1 of table 2. 8 A central constraint is imposed by the relic density selecting a thin slice in parameter space that provides a thermally averaged cross section σv FO ∼ 3×10 −26 cm 3 s −1 . We find two phenomenologically distinct regions characterised by the type of annihilation channels relevant during freeze-out: Non-resonant annihilation via Higgs exchange only leads to allowed points in the range m χ > m h /2 and m χ < m H . For points in the region m χ m H , annihilation into Higgs pairs is dominant (see the next bullet point). Points below the resonance are (mostly) excluded by the perturbativity condition (see discussion further below) as the required coupling towards small DM masses quickly becomes too large. This can be understood from eq. (3.7), representing the amplitude of the respective annihilation process. For small center-of-mass energies compared to m h /2 and m H /2, a partial cancellation takes place and the amplitude is suppressed by s. This is in contrast to the singlet scalar Higgs portal model [85][86][87] where this suppression is not present and the region below the resonance can satisfy the relic density constraint for perturbative couplings. The features at low m χ are made more apparent by zooming into a small DM mass window around the Higgs mass resonance region, m χ ∼ m h /2, as shown in figure 2.
2. Annihilation into Higgs pairs (χχ → HH, hH, hh) can be the dominant channel for m χ > m h/H , and according to the thermal momentum distribution during freeze-out, for DM masses slightly below the Higgs threshold m χ m h/H . In our scan θ ∼ 0 is preferred except for m H m h (see discussion further below). As the annihilation cross section into HH (hh) is proportional to cos 4 θ (sin 4 θ), we find that χχ → HH  figure 1. A similar behaviour was found in ref. [59]. On top of this, the observed signal strengths for a SM-like Higgs at the LHC exhibit a slight 9 As pointed out in ref. [84], the assumption of local thermal equilibrium during freeze-out can break down near the resonances. This has the effect of changing the coupling value by a factor of order O(1). However, this part of the parameter space has small v h /vs and is well beyond the sensitivity of current and future experiments. Thus, we employ the standard calculation of the DM relic density within micrOMEGAs assuming local thermal equilibrium during freeze-out. preference (around 1σ) for m H m h and sizeable θ over the SM only prediction with m h = 125 GeV. However, the presence of a second Higgs improves the fit only due to the freedom in m H . In fact, keeping m h as a free parameter as well is expected to improve the fit and broaden the 1σ CL region beyond m H m h . 10 As mentioned above, constraints from perturbative unitarity are relevant in and exclude parts of the parameter space where the measured relic density could only be matched with extremely large couplings. In the limit of small mixing, the perturbative unitarity limit in eq. (3.2) translates to JHEP04(2020)015  This limit is evident in the (m H , v h /v s )-plane for m H 10 GeV, i.e., parameter points that lie outside the boundary of the 2σ CL implies v h /v s > 713 GeV/m H . Due to the strong constraints from perturbative unitarity towards small masses, limits from the invisible Higgs decay are less relevant than e.g., in the singlet scalar Higgs portal model. We found that dropping the likelihood from invisible Higgs decay does not significantly change the results shown in figure 1.
In figure 3, we show the PLR plots for key DM observables such as the pNG DM relic density, and its annihilation cross section into SM and non-SM particles today, σv 0 . The 1σ CL region shows up as two disconnected islands. The small island at m χ m h /2 = 62.5 GeV corresponds to the h resonance, where χχ → bb channel is most dominant. The second island appears for m χ 125 GeV. As m H is profiled over, and given that m H m h is favoured, this island corresponds to the region where χχ → hh, HH is dominant -16 -

JHEP04(2020)015
during freeze-out and sets the pNG DM relic abundance to the observed value today. Note that, although σv FO ∼ 3 × 10 −26 cm 3 s −1 as required by the relic density constraint, the annihilation cross section today, σv 0 varies over many order of magnitude. This is due to the large velocity dependence of the annihilation cross section in the vicinity of a resonance and a threshold. Resonant annihilation can lead to both σv 0 smaller or larger than σv FO depending on whether the DM mass is smaller or slightly larger than the point of maximal enhancement during freeze-out. For the 1σ CL region, this behaviour can be seen in the right panel of figure 3. For annihilation into Higgs pairs, in contrast, σv 0 can only be suppressed compared to σv FO due to the smaller phase space around threshold today. Again, this behaviour can be seen for the 1σ CL region above 125 GeV in the right panel of figure 3.
The best-fit point lies in the h-resonance region exhibiting large mixing and relatively small v h /v s . The best-fit for the second Higgs mass is m H = 125.3 GeV resulting from LHC signal strength measurements. The corresponding values are summarised in table 2. Note, however, that the PLR is relatively flat within the 1σ CL region. Furthermore, as stated above, the 1σ preference for m H 125.3 GeV to some extent is a result of our choice m h = 125 GeV. This is in slight tension with the LHC Higgs signal strength and could be alleviated by treating the Higgs mass as a nuisance parameter in the fit. We therefore consider the entire 2σ CL region to be consistent with observation on a compatible level.

Marginalised posteriors
In Bayesian statistics, we rely on the Bayes' theorem: where θ are the free parameters of our model, x is the observed data, P(θ|x) is the posterior pdf, L(x |θ) is the likelihood function and π(θ) is the prior pdf. The denominator involves an integral over the free model parameters and is known as the Bayesian evidence Z. In our case of a multi-dimensional model, we are interested in 2D marginalised posterior (MP) distributions. These are constructed in the following way [88] P(θ i , θ j |x) = l = i, j dθ 1 , . . . , dθ l P(θ|x), (4.5) where we integrate over the irrelevant parameters {θ l | l = i, j}. The MP distribution above is used to define a Bayesian credible region (CR) ω in such a way that there is a probability α of containing the true values of model parameters: are less favoured as they require an extra degree of tuning of model parameters to satisfy all of the included constraints.
In comparison to the PLR plots in figure 1, the allowed regions in the MP plots are more constrained, especially where a large degree of tuning is required from marginalising over the model parameters. Again, we see a vertical stripe in the (m H , m χ )-plane. On the other hand, the second resonance region, m χ m h, H /2, is less-favoured as it falls outside the 2σ credible interval due to an extra need for tuning over v h /v s . In addition, regions where m χ > m H also appears to be fine-tuned, especially after marginalising over v h /v s .  constraint, and thus is less favoured. Lastly, in the (θ, v h /v s )-plane, the posterior mass is large for θ 0.1 rad. However, large values of θ are still allowed as they fall within the 2σ credible interval.
In figure 5, we show the MP distributions for key DM observables. In contrast to the right panel of figure 3, we do not see at least 4 orders of magnitude smaller DM annihilation cross sections than the freeze-out value; the region with velocity suppressed annihilation cross section is somehow fine-tuned and less favoured after marginalising over

Post-processing of samples
In addition to the constraints included in our global fit, we consider indirect and direct detection constraints, see sections 3.6 and 3.7, respectively. For the computation of corresponding observables, we post-process our final samples. This greatly reduces the computational time, in particular, for indirect detection constraints. The Fermi -LAT likelihood is computationally intensive due to the generation of the annihilation spectra for 2 → 2 up to 2 → 6 processes (see section 3.6). We nevertheless expect a sufficient coverage within the resulting (1 − 2)σ CL contours after combining various scans. 11 Accordingly, for the post-processed samples, we provide a frequentist interpretation only. While indirect detection constraints from Fermi -LAT observations of dSphs have a significant effect on the PLR, current direct detection experiments are not yet sensitive to our model, as we will show below. We thus refrain from including a likelihood for the latter, and restrict ourselves to comparing the model prediction to the reach of current and future experiments for this case.

Indirect detection
As explained in section 3.6, we consider two cases regarding the set of dSphs included. We take into account the likelihoods from all 45 dSphs considered in ref. [34] as well as excluding 11 As stated in footnote 7, we combine results from several MultiNest scans. This is done with various specific priors to guarantee sufficient coverage in the resonant regions as well as in regions preferred by Fermi -LAT when considering 45 dSphs. The resulting chain contains more than 3 million points.
-19 -JHEP04(2020)015 the four dSphs that show an excess, correspondingly including 41 dSphs. The latter choice only imposes an upper limit on the annihilation cross section and is described first.
In figures 6 and 7, we show the PLR plots in the planes of pNG DM model parameters as well as the DM relic abundance and annihilation cross section today, respectively, after accounting for the likelihoods of 41 dSphs. The implications on the parameter space compared to our global fit results (see figure 1) are moderate. However, for m χ 100 GeV, the Fermi -LAT limits exclude a large portion of the parameter space where the pNG DM annihilation today proceeds via χχ → HH channel, i.e., where m χ > m H . The constraint becomes stronger for smaller DM masses as lighter DM requires a larger DM number density to match the same energy density. This enhances the annihilation rate. The tendency is partly softened by the fact that for a given m H , the spectrum becomes more peaked for larger m χ , which tends to strengthen the constraints.
Taking into account the Fermi -LAT likelihood, the new best-fit point has moved to the region of dominant annihilation into Higgs pairs during freeze-out (χχ → hh in this case), see column 2 in table 2. However, m χ is slightly smaller than m h such that χχ → hh is kinematically forbidden today. Thus, σv 0 is largely suppressed as now it proceeds via (a highly off-shell) Higgs propagator (dominantly into W W * final states). Consequently, the best-fit point effectively evades any constraint from indirect detection.
Note that other indirect detection searches can impose further constraints on the parameter space. Here we would like to comment on current constraints from cosmic-ray (CR) antiproton fluxes as measured by AMS-02 [93]. While the corresponding analyses are typically plagued by large CR propagation uncertainties, recent progress has been made by fitting propagation and DM parameters at the same time [91]. This analysis provides very strong bounds on the DM annihilation cross section for DM masses above 200 GeV. While performing a respective, dedicated analysis for the considered model is beyond the scope of this work, we can, nevertheless, interpret the results of ref. [91] in parts of our model parameter space. The analysis provides limits for annihilation into a pair of Higgses with m h = 125 GeV. In our model, m H ∼ m h in the entire 1σ CL region. Moreover, for m χ > 200 GeV (where the analysis becomes constraining), the dominant annihilation channels are χχ → hh, hH and HH. Hence, for the 1σ CL region, the result from ref. [91] can be directly applied without approximation (except for neglecting the small difference between m H and m h , which is however, not expected to have a noticeable effect on the gamma-ray energy spectrum). We show the corresponding 95% CL limit as solid blue curve in the right panel of figure 7. For the 2σ CL region, in general, m H = m h . Nevertheless, the solid blue curve is expected to provide an order of magnitude estimate of the sensitivity of CR antiproton searches.
With future experiments, the pNG DM parameter space can be tested with gamma-ray observation with improved sensitivity. A large part of the region of dominant annihilation into Higgs pairs is expected to be probed in the near future by a combination of new dSphs discovered by the Large Synoptic Survey Telescope (LSST) [94] with Fermi -LAT observations. First, the inclusion of more satellite galaxies will augment the Fermi -LAT data [90,95]. Secondly, the LSST novel spectroscopic observations will provide precise measurements of J-factors, decreasing the associated astrophysical uncertainties. These improvements are expected to provide sensitivity to the thermal cross-section for DM masses up to around 600 GeV. For illustration, we show the corresponding projected limit for annihilation into bb (assuming 18 years of observation) as the dotted red curve in figure 7. Furthermore, in the right panel of figure 7 only, we display the recent projection for Galactic centre gamma-ray observations with the Cherenkov Telescope Array (CTA) [92], assuming 500 hours of exposure (no systematics), again using DM annihilation into bb as a benchmark channel, which is expected to provide a reasonable order of magnitude estimate for the sensitivity for our model. It can probe cross sections down to 5 × 10 −27 cm 3 s −1 for masses above 500 GeV, i.e., a large portion of the allowed pNG DM parameter space characterised by dominant annihilation into Higgs pairs.
In figures 8 and 9, we show the respective results after accounting for the Fermi -LAT likelihood from all 45 dSphs, i.e., including Reticulum II, Tucana III, Tucana IV and Indus II, which exhibit slight excesses with a local significance of around 2σ each [34, 69, . The solid blue curve shows the current cosmic-ray (CR) antiproton limit for the χχ → hh channel [91]. Projected limits from CTA Galactic Centre (GC) [92] (dashed green), and LSST+Fermi -LAT dSphs [90] (dotted red) for the bb channel are also shown. 70]. Interestingly, the excess can be fitted by an annihilation cross section in the ballpark of the thermal one, σv 0 ∼ 3 × 10 −26 cm 3 s −1 , i.e., in regions where the annihilation cross section today is similar to the one typically required during freeze-out. This places additional constraints on the parameter space and excludes parts of the resonant region where the cross section is highly velocity dependent, and hence σv 0 deviates strongly from σv FO . More concretely, for DM masses below the point of maximal resonant enhancement, σv 0 < σv FO and the flux today tends to be too low to fit the signal. In contrast, for DM masses above that point, the flux tends to be too high. Finally, in between, σv 0 ∼ σv FO . This is approximately the point of maximal resonant enhancement of the thermally averaged cross section during freeze-out which allows for the smallest possible couplings in the scan. If the annihilation proceeds via an on-shell Higgs (h or H), the respective cross section is proportional to Here Γ tot h/H is the total Higgs decay width, which is dominated by the partial width into SM particles if the DM mass is very close to m h/H /2, such that the corresponding phase space is suppressed. In this case, Γ tot h/H is proportional to cos 2 θ (sin 2 θ) for h (H), thereby cancelling the respective factors in the numerator of eq. (4.7). Consequently, the cross section for resonant annihilation via h or H is approximately proportional to (sin θ/v s ) 2 and (cos θ/v s ) 2 , respectively. The two regions can be recognised in the lower part of the (θ, v h /v s )-plane as the two overlapping thin bands with decreasing (∝ sin −1 θ) and increasing (∝ cos −1 θ) slope. The best-fit point falls in the second band where the pNG DM annihilation proceeds dominantly via a resonant H (see column 3 in table 2).
The second region fitting the signal is characterized by annihilation into a pair of Higgses, where (unless very close to or below threshold) no strong velocity dependence of the annihilation cross section is present, and thus σv 0 ∼ σv FO naturally. This region extends from the respective threshold up to around 200 (300) GeV within 1σ (2σ) CL region away from the best-fit point. Accordingly, the allowed region spans roughly an order of magnitude in the pNG DM mass, i.e., around 30-300 GeV. Note that in ref. [41], the pNG DM model has been considered as an explanation of the gamma-ray galactic centre excess [96][97][98][99][100][101][102][103][104] and the CR antiproton excess [105][106][107][108][109]. While a discussion of the robustness of a DM explanation of these excesses as well as an explicit interpretation is beyond the scope of this work, we briefly comment on this possibility. In particular, we distinguish two cases.
1. For pNG DM annihilation via a resonant h or H exchange in the s-channel, the composition of final states is the same as for the (singlet scalar) Higgs portal model, unless m χ > m h, H . For this case, explicit fits to the gamma-ray galactic centre excess [110] and the CR antiproton excess [107] have been performed. The latter analysis also provides a joint fit of both observations and the above-mentioned excess in the Fermi -LAT dSphs. It reveals that all three observations, if arising from DM -23 -JHEP04(2020)015 Figure 9. 2D PLR plots for the pNG DM relic abundance and its annihilation cross section today after post-processing our samples with Fermi -LAT likelihood from 45 dSphs. In the right panel, projected limits from CTA GC searches [92] and LSST+Fermi -LAT dSphs [90] are also shown.
annihilation, are compatible and point to a DM mass of around (50-60) GeV and a velocity averaged annihilation cross section today of around (1−2) × 10 −26 cm 3 s −1 .
2. The other case concerns dominant annihilation into h or H. Due to their subsequent decays into lighter SM particles, their gamma-ray spectra are typically softer. For instance, just above its threshold, the photon spectrum for χχ → HH → bbbb has the same shape as the one for χχ → bb but is shifted by a factor of 2 towards smaller energies. This is also reflected in the fact that the three observations can be fitted by DM annihilation into a SM-like Higgs for masses around the threshold [107], i.e., roughly a factor of two larger than the DM mass providing the best fit for χχ → bb.
In conclusion, the regions that fit the gamma-ray galactic centre excess and the CR antiproton excess are very similar to those preferred by the 45 dSphs (see figures 8 and 9). We expect all three observations to be well fitted by a significant subset of the parameter points fitting the 45 dSphs. Similar to the case of 41 dSphs (see figure 7), we also display the projected limits (for the bb channel) for future gamma-ray observations of dSphs from LSST+Fermi -LAT and for CTA in the right panel of figure 9. It is evident that LSST+Fermi -LAT will be able to test almost the entire 2σ CL region preferred by the 45 dSphs.

Direct detection
In figure 10, we show the PLR plots for the pNG DM-nucleon cross section at one-loop level after post-processing our MultiNest samples. These cross sections are based on the approximate expressions (left panel ) and full computations (right panel ), i.e., eqs. (3.19) and (3.20), respectively. The solid red curve shows the current sensitivity of XENON1T [4], whereas the dashed orange and dotted magenta curves show the projected sensitivities of LZ [32] and DARWIN [33], respectively.  . The solid red curve shows the current exclusion limit of XENON1T [4], whereas the dashed orange and dotted magenta curves show the projected sensitivities of LUX-ZEPLIN (LZ) [32] and DARWIN [33], respectively.
The approximate cross section overestimates the full one-loop prediction up to several orders of magnitude. For instance, parts of the 1σ CL region in the left panel are already excluded by XENON1T, while they are currently allowed when considering the full oneloop computation. In fact, we find that the entire 2σ CL region is not challenged by the current limits from XENON1T. Even the projected LZ and DARWIN experiments will probe only a small portion of the 2σ CL region. The best-fit point lies completely out of reach of these experiments. In particular, the resonance region, m χ m h /2, predicts a DM-nucleon cross section that is smaller than ∼ 10 −50 cm 2 and lies well below the proposed neutrino floor [111]. It is still interesting to see that upcoming generation of direct detection experiments are starting to probe models of DM with momentum-suppressed tree-level cross-section.

Conclusions
We performed a global fit of the pNG DM model by combining constraints from the DM relic abundance, perturbative unitarity, Higgs invisible decay, electroweak precision observables and Higgs searches at colliders. We presented our results in both frequentist and Bayesian statistical frameworks. In addition, we post-processed our samples by imposing indirect detection constraints from Fermi -LAT dwarf spheroidal galaxies within the former framework. Furthermore, we computed the one-loop pNG DM-nucleon cross sections, and compared the resulting values against the current limit from XENON1T (2018), and projected future limits from LUX-ZEPLIN (LZ) and DARWIN.
In the frequentist analysis, we found two main regions with similar profile likelihood ratio that are compatible with all observations: the Higgs funnel region where DM annihilates resonantly via one of the two Higgs bosons, m χ ∼ m h, H /2, and the region of dominant annihilation into Higgs pairs, m χ m h, H . In contrast, the region of non-resonant annihilation into SM fermions and gauge bosons is highly constrained and mostly falls outside -25 -JHEP04(2020)015 the 2σ CL region, in particular, for DM masses below the resonant region where the annihilation cross section is suppressed and requires non-perturbative couplings to match the measured relic density.
Electroweak precision observables, LEP searches and observed Higgs signal strengths at the LHC impose strong constraints on the mixing angle θ between the two Higgs bosons in our model. They require θ 0.1 rad except for the mass degenerate case m h ∼ m H , where large mixing angles are allowed as well. In fact, the observed Higgs signal strength exhibit a slight preference (around 1σ) for the latter choice. However, this preference arises from the fact that the LHC signal strengths are better fitted with a slightly heavier Higgs of around 125.3 GeV while the SM Higgs mass is fixed at 125 GeV in our scan. Hence, the fit prefers the second Higgs to have a mass of 125.3 GeV and non-suppressed couplings to the SM particles. We expect this preference to be alleviated if the SM Higgs mass is included as a nuisance parameter in the fit.
Our Bayesian results led to an even stronger constraint after marginalisation over the free model parameters. In particular, regions with a smaller volume of support fell outside the 2σ credible interval. For instance, this concerns regions where the annihilation cross section today is much larger or smaller than the canonical freeze-out cross section σv FO ∼ 3 × 10 −26 cm 3 s −1 , arising very close to the resonant condition m χ ∼ m h, H /2 or the threshold m χ ∼ m h, H . Similarly, our Bayesian results do not imply a preference for large mixing angles θ induced by Higgs signal strength observations, as it requires m H to be very close to 125.3 GeV, again, providing a small volume of support.
We computed the pNG DM-nucleon cross section at one-loop level for all of our samples after utilising the results of ref. [30]. We found that none of the points in our scan are challenged by current direct detection limits from XENON1T (2018). Future based experiments (e.g., LUX-ZEPLIN, DARWIN) will only probe a small portion of the 2σ CL region.
We took into account the Fermi -LAT likelihood by considering two different sets of dSphs. On the one hand, we considered those imposing an upper limit on the annihilation cross section only (41 dSph). The effect on the parameter space is mild and the DM mass is not constrained towards large values within the consider range. On the other hand, we considered all 45 dSphs analysed by Fermi -LAT, including the four dwarfs that show slight excesses at the level of 2σ each. These excesses can be well fitted within our model. They favour a DM mass in range 30 GeV m χ 300 GeV at the 2σ CL. We also expect a large part of this region to provide a good fit to the gamma-ray Galactic centre excess and the cosmic-ray antiproton excess seen in the AMS-02 data, if interpreted as a signal of DM annihilation.
Other indirect detection searches can further constrain our model. For instance, limits from AMS-02 antiprotons already exclude parts of the 1σ CL region in the 41-dSph fit with a DM mass around 400 GeV. Future gamma-ray observations by Fermi -LAT of newly discovered dSphs by LSST and CTA observations of the Galactic centre are expected to improve on the sensitivity and probe a significant portion of the allowed parameter space for DM masses above m h in the 41 dSphs fit. They are also expected to probe almost the entire 2σ CL region preferred by the current Fermi -LAT observations of all 45 dSphs.

JHEP04(2020)015
Acknowledgments We thank Peter Athron, Anders Kvellestad, Tim Stefaniak and Jonas Wittbrodt for helpful discussions regarding HiggsBounds/HiggsSignals. We gratefully acknowledge Da Huang for providing access to his C++ code for the calculation of one-loop pNG DM-nucleon cross section. A.S. thanks the CP3 (UCLouvain) staff for the hospitality offered during his visit where part of this work was completed.
Computational resources have been provided by the Consortium desÉquipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique

(B.5)
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.