Status of the scalar singlet dark matter model

TL;DR

This paper provides a comprehensive Bayesian and frequentist analysis of the scalar singlet dark matter model, integrating multiple experimental constraints to identify viable parameter regions and assess their naturalness.

Contribution

It offers the first combined Bayesian and frequentist global fit of the scalar singlet dark matter model using the GAMBIT framework, including extensive experimental data.

Findings

Viable parameter space for singlet masses between Higgs mass and 300 GeV.

Singlet masses above ~1 TeV can account for all dark matter.

Low-mass resonance region remains viable but is fine-tuned in Bayesian analysis.

Abstract

One of the simplest viable models for dark matter is an additional neutral scalar, stabilised by a $\mathbb{Z}_2$ symmetry. Using the GAMBIT package and combining results from four independent samplers, we present Bayesian and frequentist global fits of this model. We vary the singlet mass and coupling along with 13 nuisance parameters, including nuclear uncertainties relevant for direct detection, the local dark matter density, and selected quark masses and couplings. We include the dark matter relic density measured by Planck, direct searches with LUX, PandaX, SuperCDMS and XENON100, limits on invisible Higgs decays from the Large Hadron Collider, searches for high-energy neutrinos from dark matter annihilation in the Sun with IceCube, and searches for gamma rays from annihilation in dwarf galaxies with the Fermi-LAT. Viable solutions remain at couplings of order unity, for singlet masses between the Higgs mass and about 300 GeV, and at masses above $\sim$1 TeV. Only in the latter case can the scalar singlet constitute all of dark matter. Frequentist analysis shows that the low-mass resonance region, where the singlet is about half the mass of the Higgs, can also account for all of dark matter, and remains viable. However, Bayesian considerations show this region to be rather fine-tuned.

Figures

Tables

Table 1. Table 1 : Scalar singlet model parameters varied in our fits, along with their associated ranges and prior types.

Parameter	Minimum	Maximum	Prior
${\lambda }_{hS}$	${10}^{-4}$	10	log
$m_S$ (full-range scan)	45 GeV	10 TeV	log
$m_S$ (low-mass scan)	45 GeV	70 GeV	flat

Table 2. Table 2 : Names and ranges of Standard Model, halo and nuclear nuisance parameters that we vary simultaneously with scalar singlet parameters in our fits. We assign a flat prior to all these parameters.

Parameter		Value($\pm$Range)
Local DM density	${}^{\overline{MS}}$${\rho }_0$${}^{\overline{MS}}$	0.2–0.8 GeV cm-3
Nuclear matrix el. (strange)	${}^{\overline{MS}}$${\sigma }_s$${}^{\overline{MS}}$	$43\left(24\right)$ MeV
Nuclear matrix el. (up + down)	${}^{\overline{MS}}$${\sigma }_l$${}^{\overline{MS}}$	$58\left(27\right)$ MeV
Strong coupling	${\alpha }_s^{\overline{MS}}\left(m_Z\right)$	$0.1185\left(18\right)$
Electromagnetic coupling	$1/{\alpha }^{\overline{MS}}\left(m_Z\right)$	$127.940\left(42\right)$
Fermi coupling $\times$ ${10}^5$	${}^{\overline{MS}}$$G_{F,5}$${}^{\overline{MS}}$	$1.1663787\left(18\right)$
Higgs pole mass	${}^{\overline{MS}}$$m_h$${}^{\overline{MS}}$	124.1–127.3 GeV
Top pole mass	${}^{\overline{MS}}$$m_t$${}^{\overline{MS}}$	$173.34\left(2.28\right)$ GeV
Bottom quark mass	$m_b^{\overline{MS}}\left(m_b\right)$	$4.18\left(9\right)$ GeV
Charm quark mass	$m_c^{\overline{MS}}\left(m_c\right)$	$1.275\left(75\right)$ GeV
Strange quark mass	$m_s^{\overline{MS}}\left(2\text{GeV}\right)$	$95\left(15\right)$ MeV
Down quark mass	$m_d^{\overline{MS}}\left(2\text{GeV}\right)$	$4.80\left(96\right)$ MeV
Up quark mass	$m_u^{\overline{MS}}\left(2\text{GeV}\right)$	$2.30\left(46\right)$ MeV

Table 3. Table 3 : Parameters of each sampler for carrying out global fits of the scalar singlet model in this paper.

Scanner	Parameter	Full range	Low mass
MultiNest	nlive	20,000	20,000
	tol	${10}^{-4}$	${10}^{-5}$
Diver	NP	50,000	50,000
	convthresh	${10}^{-4}$	${10}^{-5}$
T-Walk	chain_number	$512$	$512$
	sqrtR $-$ 1	$0.01$	$0.01$
GreAT	nTrialLists	17	17
	nTrials	20,000	10,000

Table 4. Table 4 : Contributions to the log-likelihood at the best-fit point, compared to an ‘ideal’ case. The ideal is defined as the central observed value for detections, and the background-only likelihood for exclusions. Note that each likelihood is dimensionful, so its absolute value is less meaningful than any offset with respect to another point (see Sec. 8.3 of Ref.  gambit for more details of the normalisation used). The best-fit point has λ h ​ S = 6.5 × 10 − 4 subscript 𝜆 ℎ 𝑆 6.5 superscript 10 4 \lambda_{h\scriptscriptstyle S}=6.5\times 10^{-4} , m S = 62.51 subscript 𝑚 𝑆 62.51 m_{\scriptscriptstyle S}=62.51  GeV.

Likelihood contribution	Ideal	Best Fit	$\Delta \ln ℒ$
Relic density	5.989	5.989	0
LUX Run I 2015	$-$0.640	$-$0.640	0
LUX Run II 2016	$-$1.467	$-$1.468	0.001
PandaX 2016	$-$1.886	$-$1.887	0.001
SuperCDMS 2014	$-$2.248	$-$2.248	0
XENON100 2012	$-$1.693	$-$1.693	0
IceCube 79	0	0	0
$\gamma$ rays (Fermi-LAT dwarfs)	$-$33.244	$-$33.349	0.105
Higgs invisible width	0	0	0
Hadronic elements ${\sigma }_s$, ${\sigma }_l$	$-$6.115	$-$6.115	0
Local DM density ${\rho }_0$	1.142	1.142	0
$G_{\text{Fermi}}$	24.92	24.92	0
${\alpha }_{\text{EM}}$	3.350	3.350	0
${\alpha }_{\text{s}}$	6.500	6.500	0
Higgs mass	0.508	0.508	0
Top quark mass	$-$0.645	$-$0.645	0
Bottom quark mass	2.588	2.588	0
Charm quark mass	2.770	2.770	0
Light quark masses	4.844	4.844	0
Total	4.673	4.566	0.107

Table 5. Table 5 : Details of the best-fit points and posterior means, differentiated into the two main likelihood modes. Best fits are given for the case where the singlet relic density is within 1 σ 𝜎 \sigma of its observed value, and for the case where singlet particles may be a sub-dominant component of dark matter. We omit the values of the 13 nuisance parameters, as they do not deviate significantly from the central values of their associated likelihood functions.

Mode	Statistic	Relic density condition	${\lambda }_{hS}$	$m_S$ (GeV)	${\Omega }_S{h}^2$	$\log \left(ℒ\right)$	$\Delta \ln ℒ$
Low mass	Best fit	${\Omega }_S{h}^2\lesssim {\Omega }_{DM}{h}^2$	$6.5\text{×}{10}^{-4}$	$62.51$	$0.0179$	$4.566$	$0.107$
	Best fit	${\Omega }_S{h}^2\sim {\Omega }_{DM}{h}^2$	$2.9\text{×}{10}^{-4}$	$62.27$	$0.1129$	$4.431$	$0.242$
	Posterior mean	${\Omega }_S{h}^2\lesssim {\Omega }_{DM}{h}^2$	$4.3\text{×}{10}^{-3}$	$60.28$
High mass	Best fit	${\Omega }_S{h}^2\lesssim {\Omega }_{DM}{h}^2$	$9.9$	$132.5$	$1.2\text{×}{10}^{-8}$	$4.540$	$0.133$
	Best fit	${\Omega }_S{h}^2\sim {\Omega }_{DM}{h}^2$	$3.1$	$9.790\text{×}{10}^3$	$0.1131$	$4.311$	$0.362$
	Posterior mean	${\Omega }_S{h}^2\lesssim {\Omega }_{DM}{h}^2$	$3.0$	$1867$

Equations

$$\mathcal{L}=\frac{1}{2}{\mu }_S^2{S}^2+\frac{1}{2}{\lambda }_{hS}{S}^2|H{|}^2+\frac{1}{4}{\lambda }_S{S}^4+\frac{1}{2}{\partial }_{\mu }S{\partial }^{\mu }S.$$

$$m_S=\sqrt{{\mu }_S^2+\frac{1}{2}{\lambda }_{hS}{v}_0^2}.$$

$$\frac{d{n}_S}{dt}+3H{n}_S=-⟨\sigma v_{\mathrm{rel}}⟩\left(n_S^2-n_{S,\mathrm{eq}}^2\right),$$

$$⟨\sigma v_{\mathrm{rel}}⟩={\int }_{4m_S^2}^{\infty }ds\frac{s\sqrt{s-4m_S^2}{K}_1(\sqrt{s}/T)\sigma v_{\mathrm{cms}}}{16T{m}_S^4{K}_2^2(m_S/T)},$$

$$\sigma v_{\mathrm{cms}}=\frac{2{\lambda }_{hS}^2{v}_0^2}{\sqrt{s}}\frac{{\Gamma }_h(\sqrt{s})}{(s-m_h^2{)}^2+m_h^2{\Gamma }_h^2(m_h)}.$$

$${\Gamma }_{h\to SS}=\frac{{\lambda }_{hS}^2{v}_0^2}{32\pi m_h}{\left(1-4m_S^2/m_h^2\right)}^{1/2}.$$

$$N_{\mathrm{p}}=MT{\int }_0^{\infty }\varphi (E)\frac{dR}{dE}(E)dE,$$

$${\sigma }_{SI}=\frac{m_N^4}{4\pi (m_S+m_N{)}^2}\frac{{\lambda }_{hS}^2{f}_N^2}{m_h^4},$$

$$f_N=\frac{2}{9}+\frac{7}{9}\sum _{q=u,d,s}{f}_{Tq}^{(N)}.$$

$${\sigma }_l$$

$${\sigma }_s$$

$${\Phi }_i=\sum _j\frac{⟨\sigma v{⟩}_{0,j}}{8\pi m_S^2}{\int }_{E_{\text{min,i}}}^{E_{\text{max},i}}dE\frac{d{N}_{\gamma ,j}}{dE},$$

$$J_k={\int }_{\Delta {\Omega }_k}d\Omega {\int }_{\text{l.o.s.}}ds{\rho }_S^2.$$

$$\mathcal{L}(N_{\mathrm{p}}|N_{\mathrm{o}})=\frac{(b+N_{\mathrm{p}}{)}^{N_{\mathrm{o}}}{e}^{-(b+N_{\mathrm{p}})}}{N_{\mathrm{o}}!}.$$

$$\ln {\mathcal{L}}_{\text{exp}}=\sum _{k=1}^{N_{\text{dSph}}}\sum _{i=1}^{N_{\text{ebin}}}\ln {\mathcal{L}}_{ki}({\Phi }_i\cdot J_k),$$

$$\ln {\mathcal{L}}_{\text{dwarfs}}^{\text{prof.}}({\Phi }_i)=\underset{J_1\dots J_k}{\max }\left(\ln {\mathcal{L}}_{\text{exp}}+\ln {\mathcal{L}}_J\right),$$

$$\ln {\mathcal{L}}_J=\sum _{k=1}^{N_{\text{dSph}}}\ln \mathcal{N}({\log }_{10}{J}_k|{\log }_{10}{\hat{J}}_k,{\sigma }_k).$$

$${\mathcal{L}}_{{\rho }_0}=\frac{1}{\sqrt{2\pi }{\sigma }_{{\rho }_0}^{\prime }{\rho }_0}\exp \left(-\frac{\ln ({\rho }_0/{\bar{\rho }}_0{)}^2}{2{\sigma }_{{\rho }_0}^{\prime 2}}\right),$$

Full text

††journal: Eur. Phys. J. C

1. 1 Introduction

2. 2 Physics framework

3. 2.1 Model definition

4. 2.2 Relic density and Higgs invisible width

5. 2.3 Direct detection

6. 2.4 Indirect detection

7. 3 Scan details

8. 3.1 Parameters and nuisances

9. 3.2 Scanning procedure

10. 3.3 Relic density likelihood

11. 3.4 LHC Higgs likelihoods

12. 3.5 Direct detection likelihoods

13. 3.6 Indirect detection likelihoods

14. 3.7 Nuisance likelihoods

15. 4 Results

16. 4.1 Profile likelihoods

17. 4.2 Best-fit point

18. 4.3 Bayesian posteriors

19. 4.4 Vacuum stability

20. 4.5 Comparison to existing results

21. 5 Conclusion

1 Introduction

Dark matter (DM) accounts for the majority of the matter in the Universe, but its nature remains a mystery. It has been known for some time Bergstrom:2000pn ; Bertone:2004pz ; Feng:2010gw that GeV-scale particle DM can accurately reproduce the observed relic abundance of DM, provided that it has an interaction strength with standard model (SM) particles that is comparable to that of the weak force. This is the Weakly Interacting Massive Particle (WIMP) paradigm.

The simplest WIMP model is the “scalar singlet” or scalar “Higgs-portal” scenario, in which one adds to the SM a massive real scalar field $S$ uncharged under the SM gauge group SilveiraZee ; McDonald94 ; Burgess01 . $S$ is stabilised by a $\mathbb{Z}_{2}$ symmetry, and never obtains a vacuum expectation value (VEV). The only renormalisable interactions between the singlet and the SM allowed by the symmetries of the SM arise from a Lagrangian term of the form $S^{2}H^{2}$. This term gives the singlet a so-called “Higgs-portal” for interacting with the SM, leading to a range of possible phenomenological consequences. These include thermal production in the early Universe and present-day annihilation signals Yaguna09 ; Profumo2010a ; Arina11 , direct detection and $h\to SS$ decays Mambrini11 . A number of recent papers have investigated prospects for relaxing these constraints by adding additional scalars Jana16 ; Campbell16 ; Casas17 . The singlet has also been implicated in inflation Lerner09 ; Herranen15 ; Kahlhoefer15 and baryogenesis Profumo07 ; Barger09 ; JimKimmo .

The simplicity of the scenario and the discovery of the Higgs boson in 2012 Aad:2012tfa ; Chatrchyan2012 have focussed much attention on the singlet model in recent years. XENON100 and WMAP constraints were applied in Ref. Djouadi:2011aa , and an early global fit of the model using a similar range of data was performed in Ref. Cheung:2012xb . LHC Run I constraints from a CMS vector boson fusion analysis, and monojet and mono-$Z$ analyses were shown to be very weak Endo:2014cca ; indeed, monojet constraints on all minimal Higgs-portal models (i.e. scalar, fermion or vector DM interacting with the SM only via the Higgs-portal) are weak Djouadi:2012zc . Implications of the Higgs mass measurement and a detailed treatment of direct and indirect detection were explored in Ref. Cline13b , followed by the application of direct limits from the LUX and PandaX experiments Urbano:2014hda ; He:2016mls ; Escudero:2016gzx . Anti-proton data can be important in the region of the Higgs resonance Goudelis09 ; Urbano:2014hda , and competitive with the LUX limits at higher DM masses, but are ultimately prone to substantial cosmic ray propagation uncertainties. Discovery prospects at future colliders have been explored for the 14 TeV LHC and a 100 TeV hadron collider Craig:2014lda ; Han:2016gyy , and the International Linear Collider Ko:2016xwd .

The most comprehensive recent studies were presented in Refs. Cline13b ; Beniwal and Cuoco:2016jqt . The first pair of papers examined the scalar singlet scenario in light of recent (and projected) LHC Higgs invisible width measurements Belanger:2013xza ; Aad:2015pla ; Khachatryan:2016whc , the Planck relic density measurement Planck15cosmo , Planck and WMAP CMB constraints on DM annihilation at the time of recombination Cline13 ; Planck15cosmo ; Slatyer15a , Fermi-LAT analysis of gamma rays in the direction of 15 dwarf spheroidal galaxies using 6 years of Pass 8 data LATdwarfP8 , and LUX limits on the spin-independent WIMP-nucleon scattering cross-section LUX2016 . These studies also investigated the prospects for detection in gamma rays by the Cherenkov Telescope Array (CTA; Pierre14 ; Silverwood:2014yza ; Carr:2015hta ), and for direct detection by XENON1T XENON1T . Ref. Cuoco:2016jqt presented a global fit to determine the regions of the scalar singlet model space that can explain the apparent excess of gamma rays observed by Fermi towards the Galactic centre, frequently interpreted as evidence for DM annihilation Goodenough09 ; Hooper:2010mq ; Hooper:2011ti ; Abazajian:2012pn ; Gordon:2013vta ; Daylan:2014rsa ; Calore:2014nla . This included a treatment of the Planck relic density constraint, LHC invisible Higgs width constraints, direct direct search data from LUX (and projections for XENON1T and DARWIN), and constraints from Fermi-LAT searches for DM annihilation in dwarf galaxies and $\gamma$-ray lines at the Galactic centre.

Although lines and signals from the Galactic centre in the context of this model have received a reasonable amount of attention Profumo2010a ; Feng15 ; Duerr15 ; Duerr16 ; Cuoco:2016jqt , in general these signals are relevant only if the singlet is produced non-thermally, as the regions of parameter space where such signals are substantial have quite low thermal relic abundances Cline13b . Fitting the excess at the Galactic centre requires relatively large couplings, which in turn imply too little DM from thermal freeze-out. Some previous studies have solved this issue by assuming an unspecified additional production mechanism. This reduces the predictability of the theory, as the cosmological abundance of scalar singlets ceases to be a prediction. We will take a different approach, allowing for the possibility that the scalar singlet constitutes only a sub-component of DM, and permitting a different species (e.g. axions) to make up the rest. Indeed, as we show in this paper, experiments are now so sensitive to DM signals that they can probe singlet models constituting less than a hundredth of a percent of the total DM.

The purpose of the present paper is twofold. First and foremost, we provide the most comprehensive study yet of the scalar singlet scenario, in a number of ways. We augment the particle physics model parameters with a series of nuisance parameters characterising the DM halo distribution, the most important SM masses and couplings, and the nuclear matrix elements relevant for the calculation of direct search yields. These are included in the scan as free parameters, and are constrained by a series of likelihoods derived from the best current knowledge of each observable (and in some cases, their correlations). Compared to the constraints used in Refs. Cline13b ; Beniwal , we add improved direct detection likelihoods darkbit from LUX LUXrun2 , PandaX PandaX2016 , SuperCDMS SuperCDMS and XENON100 XENON2013 , as well as IceCube limits on DM annihilation to neutrinos in the core of the Sun IC79_SUSY ; IC79 . We also test some benchmark models obtained in our scan for stability of the electroweak vacuum. Given the recent preference for astrophysical explanations of the Fermi-LAT Galactic centre excess Boyarsky:2010dr ; Petrovic:2014uda ; Carlson:2014cwa ; Lee:2015fea ; Bartels:2015aea ; Cholis:2015dea ; Clark16b , we do not add this to the scan as a positive measurement of DM properties, unlike in Ref. Cuoco:2016jqt . We explore the extended parameter space in more detail than has previously been attempted, using four different scanning algorithms, and more stringent convergence criteria than previous studies. The secondary purpose of this paper is to provide an example global statistical analysis using the Global and Modular Beyond-Standard Model Inference Tool (GAMBIT) gambit , for a DM model where extensive comparison literature exists.

In Sec. 2, we describe the Lagrangian and parameters of the scalar singlet model, discuss our astrophysical assumptions, and define the nuisance parameters that we include in our global fit. Sec. 3 gives details of our scan, including the likelihood terms that we include for each constraint, the sampling algorithms we employ, and their settings. We present the latest status of the singlet model in Sec. 4, before concluding in Sec. 5.

All input files, samples and best-fit benchmarks produced for this paper are publicly accessible from Zenodo the_gambit_collaboration_2017_801511 .

2 Physics framework

2.1 Model definition

The renormalisable terms involving a new real singlet scalar $S$, permitted by the $\mathbb{Z}_{2}$, gauge and Lorentz symmetries, are

[TABLE]

From left to right, these are: the bare $S$ mass, the Higgs-portal coupling, the $S$ quartic self-coupling, and the $S$ kinetic term. Because $S$ never obtains a VEV, the model has only three free parameters: $\mu_{\scriptscriptstyle S}^{2}$, $\lambda_{h\scriptscriptstyle S}$ and $\lambda_{\scriptscriptstyle S}$. Following electroweak symmetry breaking, the portal term induces $h^{2}S^{2}$, $v_{0}hS^{2}$ and $v_{0}^{2}S^{2}$ terms, where $h$ is the physical Higgs boson and $v_{0}=246$ GeV is the VEV of the Higgs field. The additional $S^{2}$ term leads to a tree-level singlet mass

[TABLE]

Dark matter phenomenology is driven predominantly by $m_{\scriptscriptstyle S}$ and $\lambda_{h\scriptscriptstyle S}$, with viable solutions known to exist Cline13b ; Beniwal in a number of regions:

the resonance region around $m_{\scriptscriptstyle S}\sim m_{h}/2$, where couplings are very small ($\lambda_{h\scriptscriptstyle S}<10^{-2}$) but the singlet can nevertheless constitute all of the observed DM, 2. 2.

the resonant “neck” region at $m_{\scriptscriptstyle S}=m_{h}/2$, with large couplings but an extremely small relic $S$ density, and 3. 3.

a high-mass region with order unity couplings.

The parameter $\lambda_{\scriptscriptstyle S}$ remains relevant when considering DM self-interactions (e.g. Tenkanen16 ), and the stability of the electroweak vacuum. In the SM, the measured values of the Higgs and top quark masses indicate that the electroweak vacuum is not absolutely stable, but rather meta-stable Degrassi2012a . This means that although the present vacuum is not the global minimum of the scalar potential, its expected lifetime exceeds the age of the Universe. Although this is not inconsistent with the existence of the current vacuum, one appealing feature of scalar extensions of the SM is that the expected lifetime can be extended significantly, or the stability problem solved entirely, by making the current vacuum the global minimum.

The stability of the electroweak vacuum has been a consideration in many studies of scalar singlet extensions to the SM Lerner09 ; Gonderinger2010 ; Drozd2011 ; Chen2012 ; Pruna13 ; Belanger2013a ; Khan2014 ; Alanne2014 ; Han2015a ; Kanemura2015 ; Robens15 ; Robens16 , typically appearing along with constraints from perturbativity, direct detection experiments and the relic abundance of DM. As such, vacuum stability can be an interesting aspect to study of the scalar singlet model (and indeed, of any UV-complete model). In this paper however, we primarily treat the scalar singlet DM model as a low-energy effective theory, and do not consider $\lambda_{\scriptscriptstyle S}$ as a relevant parameter. In a future fit, we plan to explore renormalisation of the scalar singlet model over the full range of scales, from electroweak to Planck, including full calculations of perturbativity and the lifetime of the electroweak vacuum. Here, for the sake of interest we simply check the stability of the electroweak vacuum for a few of our highest-likelihood parameter points.

2.2 Relic density and Higgs invisible width

In order to calculate the relic density of $S$, we need to solve the Boltzmann equation Gondolo:1990dk

[TABLE]

where $n_{\scriptscriptstyle S}$ is the DM number density, $n_{{\scriptscriptstyle S},\mathrm{eq}}$ is the number density if the DM population were in chemical equilibrium with the rest of the Universe, $H$ is the Hubble rate, and $\langle\sigma v_{\mathrm{rel}}\rangle$ is the thermally averaged self-annihilation cross-section times the relative velocity of the annihilating DM particles (technically the Møller velocity). The non-averaged cross-section $\sigma v$ depends on the centre-of-mass energy of the annihilation $\sqrt{s}$, and the thermal average depends on temperature $T$. The average is given by

[TABLE]

where for convenience we have expressed the result in terms of the relative velocity of the annihilating $S$ particles in the centre-of-mass frame, $v_{\rm cms}=2\sqrt{1-4m_{\scriptscriptstyle S}^{2}/s}$. For the case of the scalar singlet model, the non-averaged cross-section for annihilation into all final states except $hh$ is Cline13b

[TABLE]

For $m_{\scriptscriptstyle S}>m_{h}$, this expression needs to be supplemented with the partial annihilation cross-section into $hh$, given in Eq. A4 of Ref. Cline13b .

We use the SM Higgs boson width $\Gamma_{h}(\sqrt{s})$ as a function of the invariant mass of the resonance $m_{h^{*}}=\sqrt{s}$, as implemented in DecayBit SDPBit .111This comes from interpolating the results contained in the tables of Ref. YellowBook13 , and does not (yet) include theoretical uncertainties or the ability to recompute the width for different values of relevant nuisance parameters, such as $\alpha_{s}$ or quark masses. Although included in DarkBit, we checked that it makes no difference to our results for $m_{\scriptscriptstyle S}<m_{h}/2$ whether or not we modify the width in the denominator of Eq. 5 corresponding to the propagator of the internal Higgs, to take into account the decay channel $h\to SS$. At tree level, the decay width of Higgs bosons to such invisible final states is

[TABLE]

This is the standard method for calculating the relic density. It assumes that kinetic decoupling of DM from other species occurs well after chemical freeze-out. If this is not the case, one must solve a coupled system of differential equations rather than the single Boltzmann equation (Eq. 3) vandenAarssen:2012ag . For the scalar singlet model, the standard approach is very accurate except at and below the Higgs resonance, where $m_{\scriptscriptstyle S}\sim m_{h}/2$. Here, the impact of a more accurate treatment on the relic density can be up to one order of magnitude in the range $53\,{\rm GeV}\lesssim m_{\chi}\lesssim 63$ GeV Binderinprep , as $\sigma v_{\mathrm{rel}}$ is resonantly enhanced, so even small values of $\lambda_{h\scriptscriptstyle S}$, where DM undergoes early kinetic decoupling, can avoid thermal overproduction. Although this effect should arguably be included for the sake of completeness in future fits, it has little impact on our final results because it only affects a relatively small mass range.

2.3 Direct detection

The predicted number of events in a direct detection experiment is

[TABLE]

with $M$ the detector mass, $T$ the exposure time, and $\phi(E)$ the detector response function. The latter encodes the fraction of recoil events of some energy $E$ that are expected to be detected, within some analysis region.

The differential recoil rate $\frac{dR}{dE}$ depends on the nuclear scattering cross-section. The scalar singlet model has no spin-dependent interactions with nuclei. The spin-independent WIMP-nucleon cross-section is

[TABLE]

where $m_{N}$ is the nucleon mass, and $f_{N}$ is the effective Higgs-nucleon coupling

[TABLE]

The three light quark nuclear matrix elements $f^{(N)}_{Tq}$ can be calculated from the nuclear matrix elements that describe the quark content of the proton and neutron,

[TABLE]

where $m_{l}\equiv(1/2)(m_{u}+m_{d})$, and $N\in\{p,n\}$. See Ref. Cline13b for details.

Halo uncertainties can have a significant impact on the interpretation of direct searches for DM Benito:2016kyp . For the DM halo in the Milky Way, we assume a generalised NFW profile, with a local Maxwell-Boltzmann speed distribution truncated at the local Galactic escape velocity. The only parameter of the density profile that we retain as a nuisance parameter is $\rho_{0}$, the local DM density, although GAMBIT makes it straightforward to also include uncertainties arising from the DM velocity profile in future fits. For the scans of this paper, we assume a most probable speed $v_{0}=235$ km s*-1* Reid:2009nj ; Bovy:2009dr , and an escape velocity of $v_{\mathrm{esc}}=550$ km s*-1* Smith:2006ym . See Ref. darkbit for details.

2.4 Indirect detection

The flux of gamma rays from DM annihilation factorises into a part $\Phi$ that only depends on the particle physics properties and a part $J$ that depends only on the astrophysical distribution of DM.

For the gamma-ray flux in energy bin $i$ with width $\Delta E_{i}\equiv E_{\text{max},i}-E_{\text{min},i}$, the particle physics factor is

[TABLE]

where $dN_{\gamma,j}/dE$ is the differential gamma-ray multiplicity for single annihilations into final state $j$, and $\langle\sigma v\rangle_{0,j}\equiv\sigma v_{j}|_{v\to 0}\equiv\sigma v_{j}|_{s\to 4m_{\scriptscriptstyle S}^{2}}$ is the zero-velocity limit of the partial annihilation cross-section into final state $j$. This is the final-state-specific equivalent of Eq. 5. We compute the partial annihilation cross-sections for the singlet model using the expressions of Appendix A of Ref. Cline13b , as implemented in DarkBit darkbit . We obtain the predicted spectra ${dN_{\gamma}}/{dE}$ for each model point by using a Monte-Carlo showering simulation, detailed in Ref. darkbit .

The astrophysics factor for a given target $k$ is

[TABLE]

Here, $\Delta\Omega_{k}$ denotes the solid angle over which the signal is integrated, l.o.s. indicates that the line element $ds$ runs along the line of sight to the target object, and $\rho_{\scriptscriptstyle S}$ is the DM mass density within it.

Neutrino telescopes also place bounds on DM models by searching for high-energy neutrinos from DM annihilation. The most likely signal in this respect comes from DM gravitationally captured by the Sun and concentrated to its core, where it would annihilate Bergstrom98b ; IC79 . This channel predominantly tests the mass and couplings of DM to nuclei rather than the annihilation cross-section, as the nuclear scattering leading to capture is the rate-limiting step for most models. Because the singlet model has no spin-dependent couplings to nuclei, neutrino telescope searches for annihilation in the Sun provide constraints only on the spin-independent scattering cross-section.

Owing to the uncertainties associated with cosmic ray propagation, we do not consider constraints from charged cosmic rays (primarily anti-protons and positrons). Radio signals coming from synchrotron emission by DM annihilation products generated in strong magnetic fields are not included in our analysis, as the associated field strengths are highly uncertain. Nor are CMB limits on DM annihilation, as Fermi dwarf limits are stronger at all masses of interest for this model. Finally, we do not consider limits implied by gamma-ray observations of the Galactic centre, whether by Fermi or ground-based gamma-ray telescopes, owing to the uncertainties involved in modelling the DM profile and astrophysical gamma-ray emission of the central Milky Way.

3 Scan details

3.1 Parameters and nuisances

A summary of the parameter ranges that we scan over for this paper is given in Tables 1 and 2.

Table 1 gives the singlet model parameters, along with the scanning priors that we use. We carry out two main types of scan: one over the full range of masses from 45 GeV to 10 TeV, intended to sample the entire parameter space, and another centred on lower masses at and below the Higgs resonance $m_{\scriptscriptstyle S}\sim m_{h}/2$, in order to obtain a more detailed picture of the resonance region.

In addition to the effect of the singlet parameters, we also consider the effects of varying a number of SM, astrophysical and nuclear parameters within their allowed experimental uncertainties. Table 2 gives the full ranges of the 13 nuisance parameters that we vary in our scans, along with their central values. We assign flat priors to all nuisance parameters in Table 2, as they are all sufficiently well constrained that their priors are effectively irrelevant.

We allow for $\pm 3\sigma$ excursions from the best estimates of the nuclear couplings. For the local DM density, we scan an asymmetric range about the central value, reflecting the log-normal likelihood that we apply to this parameter (Sec. 3.7). Detailed references for the central values and uncertainties of these parameters can be found in Ref. darkbit .

The central values of the up and down quark masses come from the 2014 edition of the PDG review PDB ; we allow these parameters to vary by $\pm 20\%$ in our fits, so as to encompass the approximate $3\sigma$ range of correlated uncertainties associated with the mass ratio likelihoods implemented in PrecisionBit SDPBit . Given the large impact that the Higgs mass can have on the phenomenology of this model, we scan an extended range for this parameter, covering more than $\pm 4\sigma$ around the central value quoted in the 2015 update to the PDG review PDG15 ($m_{h}=125.09\pm 0.24$ GeV; see Sec. 3.7). The central value and $\pm 3\sigma$ scan range for the top quark pole mass come from Ref. ATLAS:2014wva , and for all other SM nuisance parameters from Ref. PDB .

We include the local DM density and nuclear matrix elements as nuisance parameters because of their impacts on direct detection and capture of singlet particles by the Sun. The strong coupling, Higgs VEV (determined by $G_{F}$), Higgs mass and quark masses all enter into the cross-sections for annihilation and/or scattering of $S$ Cline13b . The electromagnetic coupling does not impact our fit beyond its own nuisance likelihood, but has a small effect on renormalisation of other parameters and therefore vacuum stability, which we investigate for a few benchmarks and will explore in detail in a follow-up paper.

3.2 Scanning procedure

Although 13 of the directions in the 15-dimensional parameter space are well constrained, efficiently sampling all 15 parameters simultaneously requires sophisticated scanning algorithms. We explore this space primarily with two different scanning packages interfaced via ScannerBit: a differential evolution sampler Diver, and an ensemble Markov Chain Monte Carlo (MCMC) known as T-Walk scannerbit . Both algorithms are the current state of the art when it comes to scaling with dimension scannerbit , and thus are the natural choice for this study.

Both of these algorithms are particularly well suited for multimodal distributions, and each serves a purpose in this study. T-Walk allows efficient and accurate calculation of the Bayesian posterior distribution for the target model. The package can also be used for frequentist studies if the sampling density is amplified by a judicious choice of run parameters. However, T-Walk is far less efficient at sampling the profile likelihood in high-dimensional spaces than Diver scannerbit . Because we vary 15 parameters in total, we use Diver to produce high-quality profile likelihoods. Having identified all likelihood modes, and therefore all possible locations that might meaningfully contribute to the posterior, we then use T-Walk to produce posterior distributions, checking that it does not fail to locate any of the modes identified by Diver.

In addition to the ensemble MCMC and differential evolution scans, we also combine our results with those from a more traditional MCMC, GreAT, and the nested sampling algorithm MultiNest. These are also interfaced to ScannerBit scannerbit . Although it is not typically necessary to combine results from four different algorithms, here we demonstrate the power of the GAMBIT package, which allows us to use a range of scanning procedures on the same composite likelihood, in order to produce the most robust results possible.

As discussed in Sec. 2.1, the singlet parameter space features a viable region at $m_{\scriptscriptstyle S}\approx m_{h}/2$. In this region, the annihilation of singlet DM to SM particles via s-channel Higgs exchange is resonantly enhanced, and a lower portal coupling is required to achieve the observed relic density. This region is not yet excluded by direct detection. However, probing this region of the parameter space over a large-mass range is difficult, even when using a logarithmic prior on the mass. To properly sample this region, we run a second scan with each sampler, using a flat prior over the range $m_{\scriptscriptstyle S}\in[45,70]$ GeV. We also carry out an additional specially focused low-mass scan with Diver in the “neck” region of the resonance, in order to obtain well-sampled contours in the most localised part of the allowed parameter space. We do this by excluding all points outside the range $m_{\scriptscriptstyle S}\in[61.8,63.1]$ GeV.

The convergence criteria, population size and chain details are controlled by various settings for each sampler. The settings that we use in this paper are presented in Table 3. We chose these settings after extensive testing scannerbit , to give the most stringent convergence and best exploration possible with each scanner and region. To a certain extent, some of these settings are overkill for the problem at hand, and the same physical inference could be achieved with less samples. However, the scans that we present here took only 26 000 core hours in total to compute, and the scan that dominates most of the contours (the full-range Diver scan) took just 3 hr on 10 $\times$ 24-core nodes, i.e. around 700 core hr. Compared to the time required to compute fits that include direct LHC simulations ColliderBit ; CMSSM ; MSSM , the additional sampling we do here costs practically nothing – and noticeably improves the resolution of our results. We refer the reader to Ref. scannerbit for further details of the scanners, their settings and underlying algorithms.

The profile likelihoods that we present in this paper are based on the combination of all samples from all scans, which contain $5.7\times 10^{7}$ valid samples altogether. In contrast, the posteriors that we show come exclusively from the full-range T-Walk scan.

We compute and plot profile likelihoods and posteriors using pippi pippi , obtaining profile likelihoods by maximising the log-likelihood in parameter bins over all other parameters not shown in a given plot, and posteriors by integrating the posterior density over the parameters not shown in each plot. We compute confidence regions and intervals by determining the appropriate iso-likelihood contour relative the best-fit likelihood for 1 or 2 degrees of freedom, corresponding to 1D and 2D plots, respectively. We compute Bayesian credible regions and intervals as parameter ranges containing the relevant posterior mass according to the maximum posterior density requirement. Further details can be found in Ref. pippi .

3.3 Relic density likelihood

To determine the thermal $S$ relic density for each parameter combination, we solve the Boltzmann equation (Eq. 3) numerically with DarkBit darkbit , taking the partial annihilation rates for different final states from Eq. 5 supplemented at $m_{\scriptscriptstyle S}>m_{h}$ with the expression for $\langle\sigma v\rangle_{0,hh}$ from Appendix A of Ref. Cline13b . For $m_{\scriptscriptstyle S}<150$ GeV we use the SM Higgs partial widths contained in DecayBit (from Ref. YellowBook13 ), whereas for $m_{\scriptscriptstyle S}>150$ GeV we revert to the tree-level expressions from Appendix A of Ref. Cline13b , to avoid the impact of large 1-loop corrections to the Higgs self-interaction. We determine the effective invariant rate $W_{\rm eff}$ from the partial annihilation cross-sections, and pass it on to the numerical Boltzmann solver of DarkSUSY darksusy in order to obtain $\Omega_{\scriptscriptstyle S}h^{2}$.

We implement the relic density likelihood as an upper limit only, permitting models where the thermal abundance makes $S$ a fraction of DM. Comparing with the relic abundance measured by Planck Planck15cosmo ($\Omega_{\text{DM}}h^{2}=0.1188\pm 0.0010$, at $1\sigma$), we compute a marginalised Gaussian upper limit likelihood as described in Sec 8.3.4 of Ref. gambit . Models that predict less than the measured relic density are assigned a likelihood contribution equal to that assigned to models that predict the observed value exactly. Models predicting more than the measured relic density are penalised according to a Gaussian function centred on the observed value. We adopt the DarkBit default value of 5% for the theoretical uncertainty on the relic density prediction, adding it in quadrature to the experimental uncertainty on the observed value. We note that this is a very conservative estimate of the theoretical uncertainty for this model, except in the resonance region (see Sec. 2.2).

For models that underpopulate the observed relic density, we rescale all direct and indirect signals to account for the fraction of DM that is detectable using the properties of the $S$ boson. This is internally consistent from the point of view of the model, and conservative in the sense that it suppresses direct and indirect signals in regions where the thermal abundance is less than the Planck value.

3.4 LHC Higgs likelihoods

When $m_{\scriptscriptstyle S}<m_{h}/2$, the decay $h\rightarrow SS$ is kinematically allowed, with a partial width given by Eq. 6. This is entirely invisible at hadron colliders. Constraints can be placed on the scalar singlet model parameters from measurements of Higgs production and decay rates, and the implied limit on invisible decay channels of the Higgs. For the case of SM-like couplings, the 95% confidence level upper limit on the Higgs invisible width from LHC and Tevatron data is presently at the level of 19% Belanger:2013xza . We use the DecayBit implementation of the complete invisible Higgs likelihood, based on an interpolation of Figure 8 of Belanger:2013xza .

3.5 Direct detection likelihoods

The dominant constraints on the scalar singlet model come from the LUX LUX2016 ; LUXrun2 and PandaX PandaX2016 experiments, with weaker limits also available from DarkBit based on SuperCDMS SuperCDMS and XENON100 XENON2013 . We use the DarkBit interface to DDCalc to evaluate a Poisson likelihood for observing $N_{\mathrm{o}}$ events in a given experiment, given a predicted number of signal events $N_{\mathrm{p}}$ (Eq. 7),

[TABLE]

Here $b$ is the expected number of background events in the analysis region. We model detector efficiency and acceptance effects by interpolating between values in pre-calculated tables contained in DDCalc.

3.6 Indirect detection likelihoods

The lack of evidence for anomalous gamma-ray emission from dwarf spheroidal galaxies in data collected by the Fermi-LAT experiment allows stringent constraints to be placed on the DM annihilation cross-section LATdwarfP8 . We use the Pass 8 analysis of the 6-year dataset, with the composite likelihood

[TABLE]

where $N_{\text{dSph}}$ and $N_{\text{ebin}}$ are the number of considered dSphs and the number of energy bins, respectively. The partial likelihoods $\mathcal{L}_{ki}$ are functions of the signal flux, and hence of the quantities $\Phi_{i}$ and $J_{k}$ defined in Eqs. 12 and 13, respectively.

The main results of Ref. LATdwarfP8 were obtained by profiling over the $J_{k}$ as nuisance parameters, yielding a combined profile likelihood of

[TABLE]

where

[TABLE]

Here the use of a log-normal distribution to describe the uncertainty on $J_{k}$ is a good approximation. Tabulated binned likelihoods have been provided by the Fermi-LAT experiment, and implemented in DarkBit via the gamLike package.222https://www-glast.stanford.edu/pub_data/1048/

The strongest neutrino indirect detection constraints on DM-nucleon scattering currently come from the IceCube search for annihilation in the Sun IC79 ; IC86 . We access the 79-string results via the DarkBit interface to the nulike package IC79_SUSY ; IC22Methods , which constructs a fully unbinned likelihood using event-level energy and angular information available in the published 79-string IceCube dataset, marginalised over detector systematics. We obtain predicted neutrino spectra at the Earth using WimpSim Blennow08 yield tables contained in DarkSUSY darksusy . Although IceCube limits on spin-independent scattering are not competitive with those from LUX or PandaX, for many points in the singlet parameter space they provide constraints stronger than those from SuperCDMS, and almost as strong as XENON100.

Note that the methods that we use for marginalising or profiling out additional systematic uncertainties in neutrino and $\gamma$-ray likelihoods are only applicable because the systematics are uncorrelated; the same cannot be done with a common systematic that impacts many experiments, such as the local density of DM (which affects every direct detection experiment).

The dwarf likelihood gives an identical result to what we would obtain if we were to include each of the $J$ factors as nuisance parameters in our own fit, and profile over them. The same is true of the IceCube detector systematics treated by nulike, although in that case the equivalent result would be the Bayesian one, where the corresponding nuisance parameter was marginalised over. Ideally, one would include all such nuisance parameters in the same fit, and then be free to choose at the end of a scan to profile over them all to produce profile likelihoods, or marginalise over them all to produce posterior probability densities. In practice however, the gain in accuracy achieved by doing so is generally minimal, whereas the speed gain from the ‘inline’ treatment is substantial.

3.7 Nuisance likelihoods

Following Ref. darkbit , we take the likelihood terms for the hadronic matrix elements $\sigma_{s}$ and $\sigma_{l}$ to be Gaussian, with central values and $1\sigma$ uncertainties of $43\pm 8$ and $58\pm 9$, respectively.

The canonical value of the local DM density $\rho_{0}$ is $\bar{\rho_{0}}=0.4$ GeV/cm3 (e.g. Catena:2009mf ), but this depends on assumptions such as spherical symmetry in the halo. We remain relatively agnostic with respect to this assumption by choosing a log-normal distribution for the likelihood of $\rho_{0}$, and assuming an uncertainty of $\sigma_{\rho_{0}}=0.15$ GeV cm*-3*, such that

[TABLE]

where $\sigma^{\prime}_{\rho_{0}}=\ln(1+\sigma_{\rho_{0}}/\rho_{0})$. More details can be found in Ref. darkbit .

We use the PrecisionBit implementation of SM nuisance parameter likelihoods. For the $\overline{MS}$ light quark $(u,d,s)$ masses at $\mu=2$ GeV, we use a single joint Gaussian likelihood function, combining likelihoods on $m_{u}/m_{d}$, $m_{s}/(m_{u}+m_{d})$, and $m_{s}$. We take the experimental measurements of these quantities and their uncertainties from the PDG PDB . We use Gaussian likelihoods for $G_{F}$, based on the measured value $G_{F}=(1.1663787\pm 0.0000006)\times 10^{-5}$ GeV*-2*, $\alpha_{\rm EM}$, based on the observed $\alpha_{\mathrm{EM}}(m_{Z})^{-1}=127.940\pm 0.014$ ($\overline{MS}$ scheme) PDB , and $\alpha_{s}$, using the value $\alpha_{s}(m_{Z})=0.1185\pm 0.0005$ ($\overline{MS}$ scheme), as obtained from lattice QCD. PDB . We use the quoted uncertainties as $1\,\sigma$ confidence intervals, and apply no additional theoretical uncertainty. We also apply a simple Gaussian likelihood with no theoretical uncertainty to the Higgs mass, based on the 2015 PDG result of $m_{h}=125.09\pm 0.24$ GeV PDG15 .

4 Results

4.1 Profile likelihoods

Results of our global fit analysis with all nuisances included are presented as 2D profile likelihoods in the singlet parameters in Fig. 1, and in terms of some key observables in Figs. 2 and 3. We also show the one-dimensional profile likelihoods for all parameters in red in Fig. 4.

The viable regions of the parameter space agree well with those identified in the most recent comprehensive studies Cline13b ; Beniwal . Two high-mass, high-coupling solutions exist, one strongly threatened from below by direct detection, the other mostly constrained from below by the relic density. The leading $\lambda_{h\scriptscriptstyle S}^{2}$-dependence of $\sigma_{\text{SI}}$ and $\sigma v$ approximately cancel when direct detection signals are rescaled by the predicted relic density, suggesting that the impacts of direct detection should be to simply exclude models below a given mass. However, the relic density does not scale exactly as $\lambda_{h\scriptscriptstyle S}^{-2}$, owing to its dependence on the freeze-out temperature, resulting in an extension of the sensitivity of direct detection to larger masses than might be naïvely expected, for sufficiently large values of $\lambda_{h\scriptscriptstyle S}$.333This point is discussed in further detail in Sect. 5 of Ref. Cline13b . This is the reason for the division of the large-mass solution into two sub-regions; at large coupling values, the logarithmic dependence of the relic density on $\lambda_{h\scriptscriptstyle S}$ enables LUX and PandaX to extend their reach up to singlet masses of a few hundred GeV. This is also slightly enhanced by additional $\lambda_{h\scriptscriptstyle S}^{3}$ and $\lambda_{h\scriptscriptstyle S}^{4}$ terms in $\langle\sigma v\rangle_{0,hh}$, which are responsible for the ‘kink’ seen in the border of the grey regions at $m_{\scriptscriptstyle S}\sim 600$ GeV in the left and right panels of Fig. 2.

The resonance region persists, despite being beset from all sides: invisible Higgs from above, relic density from below, indirect detection from higher masses, and direct detection from lower masses. We find a narrow “neck” of degenerate maximum likelihood directly on the resonance, with a best fit located at $m_{\scriptscriptstyle S}=62.51$ GeV, $\lambda_{h\scriptscriptstyle S}=6.5\times 10^{-4}$. The width of this region is set by a number of things:

the actual separation between the areas allowed by the invisible width and direct detection constraints, which press in from $m_{\scriptscriptstyle S}<m_{h}/2$ and $m_{\scriptscriptstyle S}>m_{h}/2$ respectively, 2. 2.

the uncertainty on the Higgs mass, which blurs the exact $m_{\scriptscriptstyle S}$ value of the resonance by $\sim$480 MeV at the level of the $2\sigma$ contours, and 3. 3.

the width of the bins into which we sort samples for plotting, which prevents anything from being resolved on scales below $\Delta m_{\scriptscriptstyle S}\sim 170$ MeV in the left panel of Fig. 1.

In addition to correctly identifying the allowed region of the parameter space, we obtain additional information from the global fit analysis beyond that seen from pure exclusion studies. Using the relic density as an upper limit, all points for which $\Omega_{\scriptscriptstyle S}h^{2}\leq\Omega_{\text{DM}}h^{2}$ have a null log-likelihood contribution, and are thus treated equally above the line in parameter space where $\Omega_{\scriptscriptstyle S}h^{2}=0.1188$. Where $\Omega_{\scriptscriptstyle S}h^{2}<\Omega_{\text{DM}}h^{2}$, we rescale the local DM density ($\rho_{0}$) as well as that in dwarfs $\rho_{\scriptscriptstyle S}$, so the direct and indirect detection likelihoods are not flat within the allowed region. Were we not to rescale signals self-consistently for the predicted relic density, the areas excluded by direct and indirect detection in the first two panels of Fig. 2 would instead closely track the standard direct and indirect sensitivity curves that many readers will be familiar with. This can be seen more clearly in Fig. 3, where we plot cross-sections rescaled by the appropriate power of $\Omega_{\scriptscriptstyle S}/\Omega_{\text{DM}}$, together with the experimental constraints from Fermi-LAT, LUX and PandaX.

Were we to instead restrict our fits to only those models that reproduce all of the DM via thermal production to with the Planck uncertainties, we would be left with a narrow band along a small number of edges of the allowed regions we have found. These edges are indicated with orange annotations in Figs. 1 and 2. At high singlet masses, the value of the late-time thermal cross-section (Eq. 4 for $T=0$) corresponding to this strip is equal to the canonical ‘thermal’ scale of 10*-26* cm3 s*-1*. At low masses, this strip runs along the lower edge of the resonance ‘triangle’ only, as indirect detection rules out models with $\Omega_{S}h^{2}=0.119$ near the vertical edge (at $m_{\scriptscriptstyle S}=62$ GeV).

In Fig. 2, we also show in grey the regions corresponding to Higgs-portal couplings above our maximum considered value, $\lambda_{h\scriptscriptstyle S}=10$, in order to give some rough idea of the area of these plots that we have not scanned (and the area that should almost certainly be excluded on perturbativity grounds were we to do so). We note that at large $m_{\scriptscriptstyle S}$, the highest-likelihood regions are all at quite large coupling values, where the annihilation cross-section is so high, and the resulting relic density is so low, that all direct and indirect signals are essentially absent – but where perturbativity of the model begins to become an issue.

4.2 Best-fit point

Our best-fit point is located within the low-mass resonance region, at $\lambda_{h\scriptscriptstyle S}=6.5\times 10^{-4}$, $m_{\scriptscriptstyle S}=62.51$ GeV. This point has a combined log-likelihood of $\log(\mathcal{L})=4.566$, shown broken down into its various likelihood components in Table 4. To put this into context, we also provide the corresponding likelihood components of a hypothetical ‘ideal’ fit, which reproduces positive measurements exactly, and has likelihood equal to the background-only value for those observables with only a limit. The overall combined ideal likelihood is $\log(\mathcal{L})=4.673$, a difference of $\Delta\ln\mathcal{L}=0.107$ with respect to our best fit. The best fit above the resonance is at $\lambda_{h\scriptscriptstyle S}=9.9$, $m_{\scriptscriptstyle S}=132.5$ GeV, with $\log(\mathcal{L})=4.540$, $\Delta\ln\mathcal{L}=0.133$.

Interpreting $\Delta\ln\mathcal{L}$ defined this way is somewhat fraught, as we do not know its distribution under the hypothesis that the best fit is correct. However, its definition is almost identical to half the “likelihood $\chi^{2}$” of Baker & Cousins Baker:1983tu , which is known to follow a $\chi^{2}$ distribution in the asymptotic limit. Our $\Delta\ln\mathcal{L}$ differs from half the likelihood $\chi^{2}$ only in that some of the components of the ideal likelihood come from the likelihood of a pure-background model, rather than from setting all predictions to their observed values. Assuming that $2\Delta\ln\mathcal{L}$ follows a $\chi^{2}$ distribution, estimating the effective number of degrees of freedom would still be difficult, as our likelihoods include many upper limits and Poisson terms, some of which have already been conditioned on the background expectation, and some of which have not. The difference between the ideal and the best-fit likelihood does nonetheless give some indication of the degree to which the Singlet DM model can simultaneously explain all data in a consistent way, and how much worse it does than the ideal model. In this sense, it gives information similar in character to the modified $p$-value method known as CLs CLs1 ; CLs2 ; Zech , which was explicitly designed for excluding models that gave poorer fits than the background model, by conditioning on the background. Were one to approximate the distribution of $2\Delta\ln\mathcal{L}$ as $\chi^{2}$ with e.g. 1–2 effective degrees of freedom, this would correspond to a rough $p$ value of between 0.6 and 0.9 in both the resonance and the high-mass region – a perfectly acceptable fit.

Next we consider parameter combinations where the singlet constitutes the entire observed relic density of DM, by restricting discussion to points with $\Omega_{S}h^{2}$ within $1\sigma$ of the Planck value $\Omega_{\text{DM}}h^{2}=0.1188\pm 0.006$ (the uncertainty includes theoretical and observational contributions added in quadrature). In this case, the best fit occurs at the bottom of the resonance, at $\lambda_{h\scriptscriptstyle S}=2.9\times 10^{-4}$, $m_{\scriptscriptstyle S}=62.27$ GeV. This point has $\log(\mathcal{L})=4.431$, which translates to $\Delta\ln\mathcal{L}=0.242$ compared the ideal model. In the high-mass region, the best fit able to reproduce the entire observed relic density is at $\lambda_{h\scriptscriptstyle S}=3.1$, $m_{\scriptscriptstyle S}=9.79$ TeV, and has $\log(\mathcal{L})=4.311$ ($\Delta\ln\mathcal{L}=0.362$). If we were to approximate the distribution of $2\Delta\ln\mathcal{L}$ as $\chi^{2}$ with 1–2 degrees of freedom, this would correspond to $p$ values of between 0.5 and 0.8 for the resonance point, and between 0.4 and 0.7 for the high-mass point. Again, these would suggest that the fit is perfectly reasonable. This indicates that there is no significant preference from data for scalar singlets to make up either all or only a fraction of the observed DM.

The four best-fit points and the corresponding relic densities are presented in Table 5.

4.3 Bayesian posteriors

By using multiple scanning algorithms in our fits, we are also able to consider marginalised posterior distributions for the singlet parameters. In Figure 4, in blue we also plot one-dimensional marginalised posteriors for all parameters, from our full-range posterior scan with the T-Walk sampler.444We choose T-Walk for this rather than MultiNest, as we find that MultiNest biases posteriors towards ellipsoidal shapes; see scannerbit for more details and example posterior maps for this same physical model. The one-dimensional posterior for $m_{\scriptscriptstyle S}$ shows that although the full-range scan has managed to detect the resonance region, this area has been heavily penalised by its small volume in the final posterior, arising from the volume effect of integrating over nuisance parameters to which points in this region are rather sensitive, such as the mass of the Higgs. The penalty is sufficiently severe that this region drops outside the $2\sigma$ credible region in the $m_{\scriptscriptstyle S}$-$\lambda_{h\scriptscriptstyle S}$ plane. We therefore focus only on the high mass modes in the righthand panel of Fig. 5, where we show the posterior from the full-range scan.

Because it is restricted to the resonance region, the low-range scan (left panel of Fig. 5) shows the expected relative posterior across this region. The fact that the resonance is so strongly disfavoured in the full-range posterior scan is an indication of its heavy fine-tuning, a property that is naturally penalised in a Bayesian analysis. This mode of the posterior accounts for less than 0.4% of the total posterior mass, indicating that it is disfavoured at almost $3\sigma$ confidence.

For the sake of understanding the prior dependence of our posteriors, we also carried out a single scan of the full parameter range with flat instead of log priors on $m_{\scriptscriptstyle S}$ and $\lambda_{h\scriptscriptstyle S}$, using MultiNest with the same full-range settings as in Table 3. Unsurprisingly, the resulting posterior is strongly driven by this (inappropriate) choice of prior, concentrating all posterior mass into the corner of the parameter space at large $\lambda_{h\scriptscriptstyle S}$ and $m_{\scriptscriptstyle S}$. The $1\sigma$ region lies above $\lambda_{h\scriptscriptstyle S}\sim 3$, $m_{\scriptscriptstyle S}\sim 3$ TeV, and the $2\sigma$ region above $\lambda_{h\scriptscriptstyle S}\sim 1$, $m_{\scriptscriptstyle S}\sim 1$ TeV.

4.4 Vacuum stability

Finally, we check vacuum stability for some interesting benchmark points.

So far, our calculations have not required any renormalisation group evolution or explicit computation of pole masses. We have simply taken the tree-level expression for $m_{\scriptscriptstyle S}$ (Eq. 2) to indicate the pole mass, and varied it and $\lambda_{h\scriptscriptstyle S}$ as free parameters. To test vacuum stability using $\overline{MS}$ renormalisation group equations (RGEs), we need to instead use these parameters along with the values of the nuisance parameters to set up boundary conditions for a set of $\overline{MS}$ RGEs. We determine values for the $\overline{MS}$ parameters that give consistent pole masses using FlexibleSUSY 555FlexibleSUSY uses numerical routines from SOFTSUSY Allanach:2001kg ; Allanach:2013kza . 1.5.1 Athron:2014yba , with SARAH 4.9.1 Staub:2008uz ; Staub:2010jh ; Staub:2012pb ; Staub:2013tta . In doing this, it becomes necessary to specify the parameter $\lambda_{\scriptscriptstyle S}$, which we set to zero at the renormalisation scale $m_{Z}$. SpecBit can then evolve the $\overline{MS}$ parameters to higher scales, using the two-loop RGEs of FlexibleSUSY, in order to test vacuum stability and also perturbativity.

For our best-fit point, the Higgs-portal coupling $\lambda_{h\scriptscriptstyle S}$ is too small to make a noticeable positive contribution to the running of the Higgs self-coupling, which reaches a minimum value of $-0.0375935$ at $2.523\times 10^{17}$ GeV. The electroweak vacuum remains meta-stable for this point, with no substantial change in phenomenology compared to the SM, where for the same Higgs and top quark masses the quartic Higgs coupling has a minimum of $-0.037631$ at $2.514\times 10^{17}$ GeV.

Next we consider a high-mass point within our $1\sigma$ allowed region: $\lambda_{h\scriptscriptstyle S}=0.5$, $m_{\scriptscriptstyle S}=1.3$ TeV. This point has a large enough coupling $\lambda_{h\scriptscriptstyle S}$ that the minimum quartic Higgs coupling is positive: $0.0522133$ at $1.40006\times 10^{9}$ GeV. We see that it is certainly possible to stabilise the electroweak vacuum within the singlet model whilst respecting all current constraints.

4.5 Comparison to existing results

The most recent study of the scalar singlet model with a $\mathbb{Z}_{2}$ symmetry and a wide range of experimental constraints was that of Beniwal et al. Beniwal . This recent study is an ideal candidate with which to compare our results, in order to check for consistency and determine the impacts of the newest experimental constraints. There are two important differences in the ingredients of our study and that of Beniwal et al. First, we include stronger DM direct detection constraints from LUX LUXrun2 and PandaX PandaX2016 , which exclude a large part of the parameter space. Second, we scan many relevant nuisance parameters, whereas previous studies have taken them as fixed. The effect of this can be seen along the boundaries of the confidence intervals, where the viable regions are always at least as large in a scan where the nuisances are allowed to vary as in one where they are fixed.

Considering these differences, we see consistency between the results of this paper and Fig. 4 of Beniwal et al. Beniwal , both in the low and high-mass parts of the $\lambda_{h\scriptscriptstyle S}$, $m_{\scriptscriptstyle S}$ parameter space. The increased size of the allowed region resulting from the variable nuisance parameters is evident along all contour edges. The behaviour of the stronger direct detection constraint is also visible, in the top left corner of the triangular part of the allowed region in the left panel of Fig. 1, and on the right side of the “neck”. In the high-mass area of the parameter space (right panel of Fig. 1), we also see LUX and PandaX cutting a large triangular region into the allowed parameter space, essentially separating the high-mass solutions into two separate likelihood modes.

5 Conclusion

The extension of the Standard Model by a scalar singlet stabilised by a $\mathbb{Z}_{2}$ symmetry is still a phenomenologically viable dark matter model, whether one demands that the singlet constitutes all of dark matter or not. However, the parameter space is being continually constrained by experimental dark matter searches. This is evident in the global fit that we have presented here, combining the latest experimental results and likelihoods to provide the most stringent constraints to date on the parameter space of this model. Direct detection experiments will fairly soon probe the entire high-mass region of the model, with XENON1T expected to access all but a very small part of each of the high-mass islands Cline13b . The resonance region will prove more difficult, though some hope certainly exists for ton-scale direct detection to improve constraints from the low-$m_{\scriptscriptstyle S}$ direction, and for future colliders focussed on precision Higgs physics to probe the edge of the region at $\lambda_{h\scriptscriptstyle S}\sim 0.02$.

We have seen that the best-fit point found in our scan does not have a notable impact on the stability of the electroweak vacuum, due to its rather small value of $\lambda_{h\scriptscriptstyle S}$. We have shown that larger values of the portal coupling can completely solve the meta-stability of the electroweak vacuum, even whilst satisfying all experimental constraints. However, the couplings required to do this are not far below where perturbativity starts to become an issue. Investigating how these competing constraints impact the allowed and preferred regions of the singlet model in a global fit will be one of the aims of our follow-up study.

In this study we have demonstrated some powerful features of the GAMBIT framework. It is now possible to easily combine likelihoods and observables from GAMBIT and existing packages in a consistent and computationally efficient way. We have varied 13 nuisance parameters in addition to the two parameters of the scalar singlet model. We have searched this parameter space using the most modern scanning algorithms available, to provide both frequentist and Bayesian statistical interpretations. By using parallel computing resources, we have achieved this with, for example, a maximum runtime for any Diver scan presented in this paper of just 3 hr. Finally, due to the modularity and flexibility of the GAMBIT system, it will be possible to include new likelihoods and/or change parts of the calculation at any time in future, in order to quickly update the analysis to take into account new experimental developments.

All input files, samples and best-fit benchmarks produced for this paper are publicly accessible from Zenodo the_gambit_collaboration_2017_801511 . The GAMBIT software is available from gambit.hepforge.org.

Acknowledgements.

We thank Jim Cline, Andrew Fowlie, Tomás Gonzalo, Julia Harz, Sebastian Hoof, Kimmo Kainulainen, Roberto Ruiz, Roberto Trotta and Sebastian Wild for useful discussions, and Lucien Boland, Sean Crosby and Goncalo Borges of the Australian Centre of Excellence for Particle Physics at the Terascale for computing assistance and resources. We warmly thank the Casa Matemáticas Oaxaca, affiliated with the Banff International Research Station, for hospitality whilst part of this work was completed, and the staff at Cyfronet, for their always helpful supercomputing support. GAMBIT has been supported by STFC (UK; ST/K00414X/1, ST/P000762/1), the Royal Society (UK; UF110191), Glasgow University (UK; Leadership Fellowship), the Research Council of Norway (FRIPRO 230546/F20), NOTUR (Norway; NN9284K), the Knut and Alice Wallenberg Foundation (Sweden; Wallenberg Academy Fellowship), the Swedish Research Council (621-2014-5772), the Australian Research Council (CE110001004, FT130100018, FT140100244, FT160100274), The University of Sydney (Australia; IRCA-G162448), PLGrid Infrastructure (Poland), Polish National Science Center (Sonata UMO-2015/17/D/ST2/03532), the Swiss National Science Foundation (PP00P2-144674), the European Commission Horizon 2020 Marie Skłodowska-Curie actions (H2020-MSCA-RISE-2015-691164), the ERA-CAN+ Twinning Program (EU & Canada), the Netherlands Organisation for Scientific Research (NWO-Vidi 680-47-532), the National Science Foundation (USA; DGE-1339067), the FRQNT (Québec) and NSERC/The Canadian Tri-Agencies Research Councils (BPDF-424460-2012).