Impact of vacuum stability, perturbativity and XENON1T on global fits of $\mathbb {Z}_2$ and $\mathbb {Z}_3$ scalar singlet dark matter

Abstract

Scalar singlet dark matter is one of the simplest and most predictive realisations of the WIMP (weakly-interacting massive particle) idea. Although the model is constrained from all directions by the latest experimental data, it still has viable regions of parameter space. Another compelling aspect of scalar singlets is their ability to stabilise the electroweak vacuum. Indeed, models of scalar dark matter are not low-energy effective theories, but can be valid all the way to the Planck scale. Using the GAMBIT framework, we present the first global fit to include both the low-energy experimental constraints and the theoretical constraints from UV physics, considering models with a scalar singlet charged under either a $\mathbb {Z}_2$ or a $\mathbb {Z}_3$ symmetry. We show that if the model is to satisfy all experimental constraints, completely stabilise the electroweak vacuum up to high scales, and also remain perturbative to those scales, one is driven to a relatively small region of parameter space. This region has a Higgs-portal coupling slightly less than 1, a dark matter mass of 1–2 TeV and a spin-independent nuclear scattering cross-section around $10^{-45}\,\hbox {cm}^2$.

1 Introduction

The discovery of the Higgs boson at the LHC makes a strong case for the existence of fundamental scalar particles. This observation immediately raises the question of whether there are other fundamental scalars that may address some of the open problems of particle physics. For example, the Standard Model (SM) of particle physics can be extended by a gauge-singlet scalar field with a stabilising symmetry in order to obtain a dark matter (DM) candidate [1,2,3]. Such a scalar singlet naturally interacts with the SM by coupling to the Higgs field and thus obtains a thermal relic abundance via the freeze-out mechanism. In spite of its simplicity, the model possesses a number of viable parameter regions that are consistent with all experimental constraints. Scalar singlets are therefore arguably the simplest realisation of the idea of weakly-interacting massive particles (WIMPs).

A remarkable feature of the scalar singlet DM model is that, just like the SM, it remains valid up to very high energies – potentially up to the Planck scale $M_\mathrm {Pl} \sim \mathcal {O}(10^{19})$ GeV. This is in sharp contrast to many alternative DM models, which are conceived only as effective low-energy theories. In fact, scalar singlets can even resolve a potential problem of the SM at high energies: for the measured values of the Higgs boson and top quark masses, the electroweak vacuum is found to be metastable, because the Higgs quartic coupling becomes negative on scales $> rsim \mathcal {O}(10^{15})$ GeV. Even though the expected lifetime of the electroweak vacuum state far exceeds the age of Universe, it is an appealing feature of scalar singlet models that the additional coupling between the Higgs and the scalar singlet affects the running of the Higgs quartic coupling at high scales and can prevent it from becoming negative [4,5,6,7,8,9,10,11,12,13,14].

In this work we present the most comprehensive study of scalar singlet DM to date by combining the information from low-energy observables, such as the relic abundance of scalar singlets and experimental constraints, with a study of the properties of the model at high energies, in particular perturbativity and vacuum stability. For this purpose we use the GAMBIT global fitting package [15], which enables the user to incorporate existing software via a backend system. Specifically, we use FlexibleSUSY 2.0.1 [16, 17] and SARAH 4.12.2 [18,19,20,21] for the renormalisation group evolution needed to study vacuum stability, DDCalc 2.0.0 [22] for DM direct detection, gamLike 1.0.0 [22] and DarkSUSY 5.1.3 [23, 24] for DM indirect detection with gamma-rays, micrOMEGAs 3.6.9.2 [25] for the relic density calculation, and Diver 1.0.4 and T-Walk 1.0.1 [26] for efficient sampling of the parameter space. This approach makes it possible to study the parameter space relevant for scalar singlets together with a number of nuisance parameters reflecting uncertainties in SM couplings and masses, the DM halo distribution, and nuclear matrix elements important for calculating nuclear scattering cross-sections. Moreover, it is easily possible to extract additional information, such as the scale at which perturbativity is violated and the expected age of the Universe in a metastable scenario.

Most studies in the past have focused on the case where a $\mathbb {Z}_2$ symmetry stabilises the singlet, and makes it a viable DM candidate [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48]. Perturbativity, vacuum stability, direct detection and the relic density of scalar singlets with a $\mathbb {Z}_3$ symmetry have also been investigated [8, 49]. The latter case introduces new phenomenology due to an additional cubic S coupling, leading to semi-annihilations. This annihilation channel can open up regions of parameter space that would otherwise be ruled out by direct detection [8], and impact indirect detection by modifying the injection spectra of light particles [50]. However, vacuum stability considerations limit the magnitude of the responsible coupling, leading to an interesting interplay of different constraints. Variants with a residual local $\mathbb {Z}_2$ or $\mathbb {Z}_3$ symmetry, arising from the breaking of a new U(1) symmetry and including an associated $Z'$ boson, have also been studied [51,52,53].

In this paper, we present an updated global analysis of the model with a global $\mathbb {Z}_2$ symmetry, and carry out the first global fit of the model with a global $\mathbb {Z}_3$ symmetry. In particular, we improve on our earlier analysis of the $\mathbb {Z}_2$ model [47] by treating the singlet self-coupling $\lambda _{\scriptscriptstyle S}$ as a free parameter, including full RGE running, considering vacuum stability and perturbativity, and incorporating the latest results from XENON1T [54, 55] and PandaX [56]. As we will see, these new direct detection results are particularly relevant, as XENON1T has sufficient sensitivity to probe the most interesting regions of parameter space. Intriguingly, rather than ruling out the model, XENON1T observes an upward fluctuation in their data, which can be interpreted as slight preference for the model that we consider.

We give details of the models in Sect. 2, of our input parameters and scanning procedure in Sect. 3, and of our observable calculations and likelihood functions in Sect. 4. The results for the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ models appear in Sects. 5 and 6, respectively. We summarise our findings in Sect. 7.

GAMBIT software can be downloaded from http://gambit.hepforge.org, and all samples, input files and best-fit points from this paper are available from Zenodo [57].

2 Model

2.1 $\mathbb {Z}_2$-symmetric model

Let us first consider the case where a real scalar singlet S is stabilised by making the Lagrangian invariant under the $\mathbb {Z}_2$ transformation $S\rightarrow -S$. The most general renormalisable scalar potential permitted by the $\mathbb {Z}_2$, Lorentz and gauge symmetries is then [58],

$$\begin{aligned} \mathrm {V}_{\mathbb {Z}_2}&= \mu _{\scriptscriptstyle H}^2 |H|^2 + \frac{1}{2} \lambda _{h}|H|^4 + \frac{1}{2}\lambda _{h\scriptscriptstyle S}S^2|H|^2 + \frac{1}{2} \mu _{\scriptscriptstyle S}^2 S^2 \nonumber \\&\quad +\, \frac{1}{4}\lambda _{\scriptscriptstyle S}S^4. \end{aligned}$$

(1)

The terms proportional to $\mu _{\scriptscriptstyle H}^2$ and $\mu _{\scriptscriptstyle S}^2$ are the Higgs and singlet bare masses, the terms proportional to $\lambda _{h}$ and $\lambda _{\scriptscriptstyle S}$ are their quartic self-couplings, and the $S^2|H|^2$ term is the portal coupling that connects the two bosons. As S never obtains a vacuum expectation value (VEV), the singlet extension is fully specified by the three parameters $\mu _{\scriptscriptstyle S}^2$, $\lambda _{h\scriptscriptstyle S}$ and $\lambda _{\scriptscriptstyle S}$. After electroweak symmetry breaking we can replace $H \rightarrow \left[ 0, (v_0+h)/\sqrt{2}\right] ^\text {T}$, with h being the SM Higgs field and $v_0 = 246$ GeV the VEV of the electroweak vacuum. Then, the portal term proportional to $\lambda _{h\scriptscriptstyle S}$ induces couplings of h to the scalar singlet S via the terms $h^2S^2$ and $v_0hS^2$. Moreover, after symmetry breaking the $\overline{MS}$ singlet mass is given by

$$\begin{aligned} \overline{m}_{\scriptscriptstyle S}= \sqrt{\mu _{\scriptscriptstyle S}^2 + \frac{1}{2}{\lambda _{h\scriptscriptstyle S}v_0^2}}, \end{aligned}$$

(2)

where $v_0$ is the $\overline{MS}$ Higgs VEV. The singlet pole mass, $m_{\scriptscriptstyle S}$, can be obtained from this using

$$\begin{aligned} m_{\scriptscriptstyle S}^2 = \overline{m}_{\scriptscriptstyle S}^2 + \mathrm {\Sigma _S}, \end{aligned}$$

(3)

where $\mathrm {\Sigma _S}$ represents loop corrections that shift the $\overline{MS}$ mass to the pole.

The only renormalisable interaction of S with the SM is through the “Higgs portal” $S^2H^2$ term. It is this term that makes it possible to have thermal production of DM in the early Universe. This portal coupling also provides potential annihilation signals [28,29,30], direct detection and $h\rightarrow SS$ decays [31]. Notice that for scalar masses less than a few TeV, the couplings $\lambda _{\scriptscriptstyle S}$ and $\lambda _{h\scriptscriptstyle S}$ necessary to explain the DM relic density remain sufficiently small to preserve perturbativity. The scalar field in this model can also feature in theories of inflation [12, 59,60,61] and baryogenesis [62,63,64].

Several viable parts of the parameter space of the scalar singlet model have yet to be probed, with the DM phenomenology essentially given by $m_{\scriptscriptstyle S}$ and $\lambda _{h\scriptscriptstyle S}$. Specifically, the parameters that have been identified to be compatible with current experimental data exist in a number of regions [34, 46, 47]:

• 1. A resonance region around $m_{\scriptscriptstyle S}\sim m_h/2$, where in spite of very small couplings ($\lambda _{h\scriptscriptstyle S}\lesssim 10^{-2}$) the singlet can nevertheless account for the entire observed relic abundancce of DM.

• 2. The resonant “neck” region at $m_{\scriptscriptstyle S}=m_h/2$, which can escape detection by the combination of large couplings and an extremely small relic S density.

• 3. A high-mass region with $\lambda _{h\scriptscriptstyle S}$ of order one.

Eventually, direct detection is expected to probe much of this remaining parameter space, leaving only large values of $\lambda _{h\scriptscriptstyle S}$ at which the theory begins to become non-perturbative [34] and a small part of the resonance region at $m_{\scriptscriptstyle S}\sim m_h/2$ untested.

Although the relic density and searches at (in)direct detection and collider experiments only probe the mass $m_{\scriptscriptstyle S}$ and the portal coupling $\lambda _{h\scriptscriptstyle S}$, the quartic self-coupling $\lambda _{\scriptscriptstyle S}$ of the scalar singlet does become relevant for the stability of the electroweak vacuum.¹ The possibility to further enlarge the expected lifetime of the electroweak vacuum or to even render it absolutely stable is an appealing feature of scalar extensions of the SM, and one of the prime motivations for our study of the scalar singlet DM model.

2.2 $\mathbb {Z}_3$-symmetric model

The symmetry group that stabilises S is not necessarily $\mathbb {Z}_2$. We will also consider a complex scalar singlet charged under a $\mathbb {Z}_3$ symmetry, with S transforming as $S\rightarrow e^{2\pi i/3}S$. This is particularly interesting because, due to the cubic $S^3$ term allowed by this symmetry, it is the simplest DM theory involving semi-annihilations [50, 66, 67], i.e. processes where two DM particles annihilate to an SM particle and another DM particle.² The most general scalar potential respecting the $\mathbb {Z}_3$ and SM symmetries is given by,

$$\begin{aligned} \mathrm {V}_{\mathbb {Z}_3}&=\mu _{\scriptscriptstyle H}^2 |H|^2 + \frac{1}{2} \lambda _{h}|H|^4 + \lambda _{h\scriptscriptstyle S}S^\dagger S|H|^2 +\mu _{\scriptscriptstyle S}^2 S^\dagger S\nonumber \\ {}&\quad +\, \lambda _{\scriptscriptstyle S}(S^\dagger S)^2 + \frac{\mu _3}{2}(S^3+S^{\dagger 3}), \end{aligned}$$

(4)

where $S^{\dagger }$ denotes the Hermitian conjugate of S. Unlike the $\mathbb {Z}_2$ model, the scalar is no longer a self-adjoint field. Instead, we have both $S^*$ and S particles, both of which contribute to the relic abundance.

This model has received significantly less attention than the $\mathbb {Z}_2$-symmetric theory, but has been studied in the context of neutrino masses [68], baryogenesis [69] and in terms of DM phenomenology [8]. The latter included constraints from vacuum stability and perturbativity along with the relic density, direct detection and invisible Higgs decays. Singlet masses below $\sim$ $53\,$GeV were ruled out by invisible Higgs decays, and the semi-annihilation process was shown to allow the model to avoid direct detection constraints in parts of parameter space where the $\mathbb {Z}_2$ model is excluded. However, as we will discuss in Sect. 4.2, vacuum stability sets a limit on $\mu _3$ and thus on the strength of semi-annihilations, so eventually this model also comes within reach of tonne-scale direct detection experiments.

3 Input parameters and sampling

3.1 Parameters and nuisances

In Ref. [47], we studied the direct phenomenological implications of the $\mathbb {Z}_2$ symmetric scalar singlet model defined at a low energy scale, without considering renormalisation of the theory, running couplings nor vacuum stability. In this sense, Ref. [47] treated the scalar singlet as an effective field theory at the scale of the scalar mass. In this study, we will go on to examine the implications of considering the scalar singlet as a UV-complete theory. The input parameters and their required ranges are necessarily different for each of these studies.

The parameters and ranges that we scan over in our fits, along with those that we hold fixed, are presented in Tables 1–3.

Table 1 Model parameters that we vary in our fits, as well as the ranges over which we vary them, and the types of priors that we apply to the sampling. The mixed prior for the parameter $\mu _3$ consists of two separate scans. One scan employs a flat prior between 0 and 1 GeV and a logarithmic prior from 1 to 4 TeV, whereas the other scan employs a flat prior for the full range

Table 1 gives the parameters of the scalar singlet models and the priors on them that we adopt in our scans. We carry out two main types of scans: the first considering masses across the entire parameter space, from 45 to 10 TeV, and a second focussed on masses at and below the Higgs resonance $m_{\scriptscriptstyle S}\sim m_h/2$, in order to obtain better sampling of this region. Notice that we do not scan over DM masses below 45 GeV, as this part of parameter space is robustly excluded by the combination of direct detection searches, constraints on the invisible decay width of the Higgs, and the singlet relic density. Note also that although the $\overline{MS}$ mass $\overline{m}_{\scriptscriptstyle S}$ (at scale $\overline{m}_{\scriptscriptstyle S}$) is the actual input parameter in our scans, our effective prior range is defined in terms of the S pole mass, as we scan over a larger range of $\overline{m}_{\scriptscriptstyle S}$ but apply a cut on the S pole mass after spectrum generation.

To study the phenomenology of a given model, one must be able to compute perturbative expressions, such as pole masses and loop-corrected scattering cross-sections. We thus demand perturbativity as part of the likelihood analysis, by invalidating points in parameter space where any of the dimensionless coupling parameters exceed $\sqrt{4\pi }$. The choice of this value is similar to that in other studies [59, 70].

Table 2 Names and ranges of SM, nuclear and halo nuisance parameters that we vary simultaneously with scalar singlet parameters in our fits. We sample all these parameters using flat priors

In addition to the scalar singlet parameters, we also vary a number of nuclear, SM and astrophysical parameters within their allowed experimental or observational uncertainties. Table 2 gives the full ranges of all the nuisance parameters that we consider, along with the central values that we adopt. We use flat priors for sampling all nuisance parameters, as each of the parameters is well enough constrained that the choice of prior has no effect.

In this paper, where we include renormalisation of the input masses, we trade the Higgs pole mass for the $\overline{MS}$ mass $\overline{m}_h= \sqrt{-2\mu _{\scriptscriptstyle H}^2}$, defined at the scale $\overline{m}_{\scriptscriptstyle S}$. We then compute the physical pole mass from the input parameters (see discussion in Sect. 4.1). This relationship is affected by radiative corrections from the scalar singlet mass, so the relationship between $\overline{m}_h$ and the pole mass is not constant throughout the parameter space, and we must therefore scan a large range for $\overline{m}_h$. The resultant value for the pole mass $m_h$ is constrained by the likelihood function described in Sect. 4.7.

Table 3 Names and values of parameters that we hold fixed in our fits

We scan over a range of $\pm 3\sigma$ around the best estimates of the strong coupling, top pole mass, nuclear matrix elements, the most probable DM speed in the Milky Way halo $v_\text {mean}$, and the Galactic escape velocity at the solar position $v_\text {esc}$. We use the same parameter to control $v_\text {mean}$ and the rotation speed, $v_{\text {rot}}$ of the galactic disk, as these can be taken as approximately equal under the assumption of a smooth, spherical DM halo. We apply a log-normal likelihood to the local DM density $\rho _0$, so we scan an asymmetric range about the central value for this parameter. Details of the likelihoods that we apply to these parameters, along with references for their central values and measured uncertainties, can be found in Sect. 4.7.

We scan over the nuclear matrix elements and local DM density because they each have a significant impact on direct detection. The strong coupling and Higgs mass enter into the cross-sections for annihilation and nuclear scattering of S [34]. In Ref. [47], we included 13 nuisance parameters. In that study, we determined that varying the masses of the bottom, charm, strange, up and down quarks, the Fermi coupling and the electromagnetic couplings within their experimentally-allowed ranges did not have any significant effect on the results. Here we therefore fix those parameters (Table 3). On the other hand, including the uncertainties of the local DM velocity profile would have a more important effect. We therefore also include the most probable DM speed $v_\text {mean}$ and the local Galactic escape speed $v_\mathrm {esc}$ as nuisance parameters in our fits here. This results in a total of 8 nuisance parameters, or 11 and 12 parameters in total for our respective scans of the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ models.

The reduction in the total number of nuisance parameters here compared to Ref. [47] is also intended to counter-act the increased computational requirements for this global fit. The likelihood is significantly more demanding of computing resources due to the need to solve the RGEs and compute pole masses, and as a result takes longer to compute. We have also replaced the relatively small prior on the Higgs pole mass in Ref. [47] with a much less constrained $\overline{MS}$ mass, in order to be able to effectively sample Higgs masses around the observed value across the whole scalar singlet parameter space. Therefore, although we have fewer nuisance parameters in these global fits, they actually require more computational resources than those of Ref. [47].

Our adopted masses for the Z boson, $\tau$ lepton and the other quarks, as well as our chosen Fermi and electromagnetic couplings, come from the 2017 compilation of the Particle Data Group [71] (Table 3).

3.2 Scanning procedure

Although many directions in parameter space are well constrained, efficient sampling of both the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ scalar singlet models still requires sophisticated sampling algorithms. We scan the parameter space with a differential evolution sampler Diver [26], and an ensemble Markov Chain Monte Carlo (MCMC) T-Walk [26]. Both algorithms are particularly well-suited to multi-modal problems in many dimensions. Diver is an optimiser best suited to mapping the profile likelihood, whereas T-Walk is better suited to obtaining the Bayesian posterior. Using Diver, we first obtain well-sampled profile likelihoods, and make sure to identify all modes of the likelihood surface. This provides information about the locations able to potentially contribute to the posterior. We then obtain posterior distributions using T-Walk, making sure that it does not fail to identify any of the modes found by Diver.

Sampling the resonance region $m_{\scriptscriptstyle S}\approx m_h/2$ can be challenging when scanning over a large mass range. We run an additional Diver scan in this region to ensure sufficient sampling, employing a flat prior between $m_{\scriptscriptstyle S}=45$ and 70 GeV. The “neck” part of the resonance is even more difficult; here we perform a third, even more focussed scan, excluding any points with S pole masses not in the range $m_{\scriptscriptstyle S}\in [61.8,63.1]$ GeV.

We repeat all scans of the $\mathbb {Z}_3$ model with both flat and log priors on $\mu _3$, as the choice of prior on this parameter can have a significant impact on the completeness with which the profile likelihood of this model is sampled.

We perform identical scans with and without the requirement that the singlet fully stabilises the electroweak vacuum. With this requirement imposed, models with a metastable electroweak vacuum (such as the SM) are excluded. Although such models are not physically invalid, the ability to make the electroweak vacuum absolutely stable is theoretically appealing. We do not perform a scan over the low-mass range with this additional constraint, as it is ruled out except for the very top of the neck region ($\lambda _{h\scriptscriptstyle S} > rsim 0.2$). We also carry out additional scans with Diver using the older 2017 XENON1T constraint [54] instead of the recent 2018 result [55], for comparison.

The population and convergence settings that we use for each sampler are given in Table 4. These settings are based on a series of extensive tests and optimisations [26]. The Diver scans that we present here each used 3400 Intel Xeon Phi 7250 (Knights Landing) cores, for approximately 86 h in total across all scans. For the 6 T-Walk scans, instead of using a fixed tolerance associated with the parameter, we found that more reliable sampling could be obtained in the current study by simply running on 1360 cores and halting scans after 23 h, using the parameter newly implemented in T-Walk 1.0.1.

The posteriors that we show come from the T-Walk scans only. Our profile likelihood plots are based on the final merged set of samples from all scans with common physical requirements and likelihoods. This includes both Diver scans and any relevant T-Walk scans, any targeted low-mass or neck scans, and scans with different priors on $\mu _3$. Without the requirement of absolute vacuum stability, our final profile likelihoods (i.e. including XENON1T 2018 results) are based on a total of $4.9\times 10^7$ and $1.9\times 10^8$ samples for the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ models, respectively. With the requirement of absolute vacuum stability, the profile likelihoods are based on $3.3\times 10^7$ and $6.4\times 10^7$ samples for the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ models, respectively.

We produce posteriors and profile likelihoods with pippi [72], basing our posteriors on the maximum posterior density requirement.

Table 4 Sampling parameters for global fits of the $\mathbb {Z}_2$- and $\mathbb {Z}_3$-symmetric scalar singlet models in this paper

4 Physics framework and likelihood details

4.1 Pole masses and $\overline{MS}$ parameters

To investigate vacuum stability and perturbativity, and to calculate observables contributing to the likelihood, we require pole masses and running parameters consistent with known SM data. We obtain these using the two-loop RGEs of FlexibleSUSY 2.0.1 [16, 17], via the SpecBit [73] interface within GAMBIT [15]. FlexibleSUSY uses SARAH [18,19,20,21], along with parts of SOFTSUSY [74, 75] and other higher-order corrections [76,77,78,79,80,81,82,83,84].

As we vary the scalar singlet mass over two orders of magnitude, we use the EFTHiggs mode of FlexibleSUSY, which uses the algorithm developed in Ref. [76] and refined in Ref. [17]. This implements a matching and running procedure for effective field theories, which is appropriate when $\overline{m}_{\scriptscriptstyle S}\gg m_t$, while not compromising the precision of the Higgs pole mass calculation (due to normal EFT uncertainties from missing $\mathcal {O}(p^2 /\overline{m}_{\scriptscriptstyle S}^2)$ terms) when $\overline{m}_{\scriptscriptstyle S}$ is close to $m_t$.

We use two-loop SM RGEs between the electroweak scale and the scale of new physics. We take the scale of new physics to be the scalar singlet running mass $\overline{m}_{\scriptscriptstyle S}$. At $\overline{m}_{\scriptscriptstyle S}$ we perform a matching between the SM and the scalar singlet DM model³ using the FlexibleEFTHiggs matching conditions given in Ref. [17]. At scales larger than $\overline{m}_{\scriptscriptstyle S}$, we use two-loop RGEs for the scalar singlet model.

For the $\mathbb {Z}_2$ model, the inputs to the FlexibleSUSY spectrum generator are the $\overline{MS}$ Lagrangian parameters $\lambda _{h\scriptscriptstyle S}$, $\lambda _{\scriptscriptstyle S}$, $\mu _{\scriptscriptstyle H}$ and $\mu _{\scriptscriptstyle S}$, defined at the renormalisation scale $Q=\overline{m}_{\scriptscriptstyle S}$. For the $\mathbb {Z}_3$-invariant version, this parameter set is extended to include $\mu _3(\overline{m}_{\scriptscriptstyle S})$. The parameters $\lambda _{h\scriptscriptstyle S}(\overline{m}_{\scriptscriptstyle S})$, $\lambda _{\scriptscriptstyle S}(\overline{m}_{\scriptscriptstyle S})$ and $\mu _3(\overline{m}_{\scriptscriptstyle S})$ are obtained directly from the model parameters sampled by ScannerBit, while $\mu _{\scriptscriptstyle H}(\overline{m}_{\scriptscriptstyle S})$ and $\mu _{\scriptscriptstyle S}(\overline{m}_{\scriptscriptstyle S})$ are obtained by inverting $\overline{m}_h^2 = -2\mu _{\scriptscriptstyle H}^2$ and Eq. 2 respectively.⁴ FlexibleSUSY fixes the remaining $\overline{MS}$ parameter $\lambda _{h}$ to ensure correct electroweak symmetry breaking. It also accepts additional input data in the form of an SMINPUTS block, as defined in the second SUSY Les Houches Accord (SLHA2) [85]. These are given in Tables 2 and 3.

Once FlexibleSUSY has determined a consistent set of $\overline{MS}$ parameters, it computes the singlet and Higgs pole masses, using

$$\begin{aligned} m_h^2 = -2\mu _{\scriptscriptstyle H}^2 + \mathrm {\Sigma _H}, \end{aligned}$$

(5)

and Eq. 3. Here $\mathrm {\Sigma _S}$ and $\mathrm {\Sigma _H}$ incorporate the full one-loop self-energies generated with the help of SARAH 4.12.2 [18,19,20,21]. $\mathrm {\Sigma _H}$ additionally includes all two-loop contributions from dominant orders, $\mathcal {O}(\alpha _t^2 + \alpha _t \alpha _s)$. Note that in this approach, we obtain the Higgs pole mass as an output rather than an input parameter, and scan the parameter space by varying the input $\overline{MS}$ mass $\overline{m}_h$. As the value of the scalar singlet mass can have a significant impact on the relationship between $\overline{m}_h(\overline{m}_{\scriptscriptstyle S})$ and $m_h$, we allow $\overline{m}_h$ to vary from 80–180 GeV. This is sufficient to permit a value of $\overline{m}_h$ that can give a 125 GeV pole mass throughout the scalar singlet parameter space. We penalise all other points using a Gaussian likelihood centred on the experimentally-measured mass ($m_h= 125.09 \pm 0.24$ GeV; see Sect. 4.7).

Note that in this setup neither the Higgs pole mass $m_h$, nor the quartic coupling $\lambda _{h}$ are inputs. Nonetheless, $m_h$ is well constrained by the likelihood, so it is important that we can calculate both $m_h$ and $\lambda _{h}$ consistently using the procedure described above.

4.2 Vacuum stability and perturbativity

With the running $\overline{MS}$ parameters of the scalar singlet model obtained as described in the previous section, it is now possible to run the couplings from the electroweak scale to the Planck scale, $M_{\text {Pl}} = 1.22\times 10^{19}$ GeV and test for vacuum stability.

We classify the stability of the electroweak vacuum in three possible ways:

Stable :

  If $\lambda _{h}(Q)>0$ for all $Q<M_{\text {Pl}}$, the electroweak vacuum is the global minimum of the Higgs potential for all Q up to the Planck scale, and is therefore absolutely stable with respect to quantum fluctuations.

Metastable :

  If $\lambda (Q_0)<0$ for any $Q_0<M_{\text {Pl}}$, but the electroweak vacuum has an expected lifetime that exceeds the age of the Universe.

Unstable :

  If $\lambda (Q_0)<0$ for any $Q_0<M_{\text {Pl}}$ and the electroweak vacuum has an expected lifetime that is less than the age of the Universe.

To distinguish between the latter two cases, and incorporate a likelihood penalty associated with vacuum decay, one should calculate the decay rate of the electroweak vacuum. A detailed description of how to obtain the decay rate and estimate the probability that the vacuum would decay within the age of the Universe can be found in section 2.5 of Ref. [73] and references therein.

At large field values, the potential can be approximated as $V \approx \frac{1}{4} \lambda _{h}|H|^4$. This can be used to find the so-called “bounce” solution to the Euclidean equation of motion, and obtain the bounce action [86],

$$\begin{aligned} B=\frac{8\pi ^2}{3\left| \lambda _{h}\right| }. \end{aligned}$$

(6)

The rate of bubble nucleation per unit volume per unit time can be estimated from the bounce action using,

$$\begin{aligned} \varGamma \approx \varLambda _B^4 e^{-B}, \end{aligned}$$

(7)

where $\varLambda _B$, the scale at which $\lambda _{h}$ is minimised, has been introduced following Ref. [87]. As we are interested in the probability that the Universe would have decayed in our past light cone, we introduce the lifetime of the Universe, $T_U \approx e^{140} / M_{\text {Pl}}$, and use it to define the volume of the past lightcone, $T_U^4$. The predicted number of decays in our past lightcone is therefore $T_U^4\varLambda _B^4e^{-B}$. We can hence define a likelihood contribution for no decay having occurred in our past light cone as,

$$\begin{aligned} \mathcal {L}=\exp \left[ -\left( e^{140}\frac{\varLambda _B}{M_{\text {Pl}}}\right) ^4\exp \left( {-\frac{8\pi ^2}{3\left| \lambda _{h}(\varLambda _B)\right| }}\right) \right] , \end{aligned}$$

(8)

based on the Poisson probability of the Universe having decayed out of the electroweak vacuum by the present day.

The actual predicted lifetime in years is,

$$\begin{aligned} \frac{\tau }{\text {yr}} = 2.09\times 10^{-32} \left( \frac{\text {GeV}}{r}\right) , \end{aligned}$$

(9)

where the rate r is given by,

$$\begin{aligned} r \approx T_U^3\varLambda _B^4 e^{-B}. \end{aligned}$$

(10)

For the SM, this gives a predicted lifetime of $\sim 1.1\times 10^{99}$ years.

Because the dominant contribution from a scalar singlet to the running of $\lambda _{h}$ is always positive, the electroweak vacuum can only become more stable in the models we consider in this paper than it is in the SM. As the probability of vacuum decay is already very small even in the metastable SM, the effect of going from a metastable vacuum to an absolutely stable one has a negligible impact on the composite likelihood.⁵ However, because the scenario of absolute stability is theoretically appealing, we repeat our global fits with the strict condition that all models must be absolutely stable, invalidating all parameter combinations that give a metastable vacuum.

In the $\mathbb {Z}_3$ model, low-scale vacuum stability gives an additional constraint on the $\mu _3$ parameter. If $\mu _3$ is large, the scalar potential can posess $\mathbb {Z}_3$-breaking minima already at the weak scale, which would be degenerate with or deeper than the SM vacuum. This can be avoided by placing an upper bound on the $\mu _3$ parameter. We adopt the condition given in Ref. [8] for an absolutely stable SM vacuum, as an upper limit on $\mu _3$:

$$\begin{aligned} \mu _3 \le 2\sqrt{\lambda _{\scriptscriptstyle S}}m_{\scriptscriptstyle S}\, . \end{aligned}$$

(11)

This constraint can be relaxed slightly by allowing for the possibility of a $\mathbb {Z}_3$-breaking minimum with a lower potential energy than the SM vacuum, but an SM vacuum with a decay half-life longer than the age of the Universe (see Ref. [8]). We do not consider this possibility, as part of our interest in studying scalar singlet DM, particularly in this global fit, is the appeal of removing metastability from the SM altogether.

We will also require that the scalar singlet couplings remain positive, such that the scalar singlet potential is bounded from below. This means that we can isolate our study of the electroweak vacuum to the Higgs dimension only. This analysis neglects the possibility of a second minimum forming in the S direction of the potential, which is possible when $\mu _S^2<0$ and $\lambda _{h\scriptscriptstyle S}$ is sufficiently large [29]. Due to the nature of the RGEs for the dimensionless scalar couplings, $\lambda _{h\scriptscriptstyle S}$ and $\lambda _{\scriptscriptstyle S}$, these couplings only grow with scale.

We let $\varLambda _P$ denote the scale where the dimensionless couplings become larger than our upper bound for perturbativity $\sqrt{4\pi }\approx 3.54$. If $\varLambda _P<m_{\scriptscriptstyle S}$, we invalidate the point; otherwise, we record the scale $\varLambda _P$ for later analysis.

There is an important caveat to our definition of vacuum stability and how we apply this as a constraint on the parameter space. In many cases, increasing the values of the dimensionless couplings in the scalar singlet sector ($\lambda _{h\scriptscriptstyle S}$ and $\lambda _{\scriptscriptstyle S}$) results in the theory becoming non-perturbative at energy scales as low as the electroweak scale. Because perturbation theory is no longer applicable in this case, we cannot compute the running of the quartic Higgs coupling to the typical scales of instability, so our analysis does not encounter a minimum and thus renders the electroweak vacuum “stable”. Such parameter combinations therefore pass the test for stability. This caveat is acceptable, because such models can still be filtered out (if desired) based on the extremely low scale at which perturbativity is broken, as given by $\varLambda _P$. Nevertheless, it is important to consider the order in which we apply these constraints when interpreting the results in Sects. 5 and 6.

Fig. 1

The diagrams for annihilation, semi-annihilation, scalar-nucleon scattering and Higgs invisible decays in the $\mathbb {Z}_3$ scalar singlet model. Here N denotes nucleons, f are SM fermions and V SM gauge bosons. Except for the semi-annihilation processes, the equivalent diagrams apply in the $\mathbb {Z}_2$ scenario but with $S^*$ replaced with S

4.3 Relic density

In the early Universe, scalar singlet DM would have been in thermal equilibrium with SM particles. The annihilation processes in the top two rows of Fig. 1 would have occurred frequently compared to the Hubble expansion rate. As the Universe expanded and cooled, the density of the scalar fields would have decreased, making forward annihilation reactions extremely rare and causing the DM abundance to freeze out. To find this abundance, we solve the Boltzmann equation [88]

$$\begin{aligned} \frac{dn_{\scriptscriptstyle S}}{dt}+3Hn_{\scriptscriptstyle S}=-\langle \sigma v_\mathrm {rel}\rangle \left( n_{\scriptscriptstyle S}^2-n_{{\scriptscriptstyle S},\mathrm {eq}}^2 \right) \,. \end{aligned}$$

(12)

Here H is the Hubble rate, $\langle \sigma v_\mathrm {rel}\rangle$ is the thermal average of the relative velocity of DM particles times their self-annihilation cross-section, the number density of DM is given by $n_{\scriptscriptstyle S}$, and its equilibrium number density by $n_{{\scriptscriptstyle S},\mathrm {eq}}$.

When semi-annihilation processes are possible, as in the $\mathbb {Z}_3$ scalar singlet model, then Eq. (12) must be modified. The tree-level semi-annihilation processes, where two DM particles can annihilate to a DM particle and an SM particle, are shown in Fig. 1. In the $\mathbb {Z}_3$ model, the relic abundance consists of equal parts $S^*$ and S, as each annihilation processes requires both an S and $S^*$, and the semi-annihilation process can occur equally rapidly via $SS\rightarrow S^*h$ and $S^*S^*\rightarrow Sh$. We can therefore treat S and $S^*$ as the same particle in the Boltzmann calculation, by including a factor of 1 / 2 [8].

$$\begin{aligned} \begin{aligned} \frac{dn_{\scriptscriptstyle S}}{dt}+3Hn_{\scriptscriptstyle S}=&-\langle \sigma v_\mathrm {rel}\rangle \left( n_{\scriptscriptstyle S}^2-n_{{\scriptscriptstyle S},\mathrm {eq}}^2 \right) \\ {}&-\frac{1}{2}\langle \sigma v_\mathrm {rel}\rangle _{SS\rightarrow hS}\left( n_{\scriptscriptstyle S}^2-n_{\scriptscriptstyle S}n_{{\scriptscriptstyle S},\mathrm {eq}} \right) , \end{aligned} \end{aligned}$$

(13)

where $\langle \sigma v_\mathrm {rel}\rangle$ is the thermally averaged self-annihilation cross-section without semi-annihilations, and $\langle \sigma v_\mathrm {rel}\rangle _{SS\rightarrow hS}$ is the equivalent for the semi-annihilation channel. We define a semi-annihilation fraction

$$\begin{aligned} \alpha =\frac{1}{2}\frac{\langle \sigma v_\mathrm {rel}\rangle _{SS\rightarrow hS}}{\langle \sigma v_\mathrm {rel}\rangle +\frac{1}{2}\langle \sigma v_\mathrm {rel}\rangle _{SS\rightarrow hS}}, \end{aligned}$$

(14)

which we record for each sampled point in the $\mathbb {Z}_3$ parameter space. To deal with semi-annihilations in the $\mathbb {Z}_3$ model, we compute $\varOmega _{\scriptscriptstyle S} h^2$ using micrOMEGAs 3.6.9.2 [25], with the setting .

As in Ref. [47], we employ the measured relic density $\varOmega _\text {DM} h^2 = 0.1188\pm 0.0010$ [89] as an upper limit, allowing models where the S relic abundance indicates that it is only a fraction of the observed DM. We use a marginalised Gaussian upper limit likelihood [15] for this purpose, adopting the default 5% theoretical uncertainty offered by DarkBit and combining it in quadrature with the uncertainty on the measured value. We self-consistently rescale all direct and indirect signals for the thermal S relic density at each point in the parameter space.

4.4 Direct detection

Scalar singlet DM is strongly constrained by limits on the DM-nucleon scattering cross-section from direct detection. The corresponding tree-level processes are represented in the bottom left diagrams of Fig. 1.

We apply direct detection constraints using DarkBit, drawing on the DDCalc [22] implementations of the experimental results of LUX [90], PandaX [56, 91] and XENON1T [54, 55]. We emphasise that since all three experiments have similar sensitivity, a consistent combination of the respective likelihoods is essential to infer accurate constraints on the parameter space.

For a given experiment, the likelihood of observing N direct detection events, given a predicted number of signal events $N_{\mathrm {p}}$, follows a Poisson distribution

$$\begin{aligned} \mathcal {L}(N|N_{\mathrm {p}}) = \frac{(b+N_{\mathrm {p}})^{N} \, e^{-(b+N_{\mathrm {p}})}}{N!} \; , \end{aligned}$$

(15)

where b denotes the expected number of background events within the analysis region. We interpolate between values in pre-calculated tables contained in DDCalc in order to determine the detector efficiencies and acceptance effects. The likelihood in Eq. (15) is then obtained by recasting the experimental results contained in DDCalc [22] for each experiment.

In particular, for the recent XENON1T results [55] we employ the data collected within the core mass of 0.65 t, which in comparison to the full 1.3 t dataset has a substantially smaller level of surface and neutron background events. We take the total detection efficiency as a function of the recoil energy from Ref. [55], weighting it by an additional factor 0.65 / 1.3. Moreover, we only consider events within the reference region defined as the area between the median of the nuclear recoil band and the $2\sigma$ quantile, leading to an additional factor 0.475 in the detection efficiency. In our analysis, we then divide the events into two energy bins, based on a separation of the S1 signal into the intervals $[3\,\text {PE},\,35\,\text {PE}]$ and $[35\,\text {PE},\,70\,\text {PE}]$. To this end, we convert a given nuclear recoil energy into an expected S1 signal using Fig. 3 of Ref. [55], and determine its probability to fall in either of the S1 bins by assuming S1 to be Poisson distributed. Furthermore, we assume that the electron recoil background is constant in S1, while we take the energy dependence of the neutron background from Ref. [92]. Assuming for simplicity that the remaining (sub-dominant) background contributions fall into the first energy bin, this gives 0.46 and 0.34 expected background events for the lower and upper energy bins of our analysis, respectively, compared to 0 and 2 observed events. With these assumptions, we obtain a 90% CL upper bound on the spin-independent scattering cross-section of DM in good agreement with the published bound, and also reproduce the slight preference (less than $2 \sigma$) for a non-zero cross-section at large DM masses.

4.5 Indirect detection

Searches for anomalous gamma-ray emission in dwarf spheroidal galaxies constrain the DM annihilation cross-section. The expected flux of gamma rays is,

$$\begin{aligned} \varPhi _i = \sum _j\frac{\langle \sigma v\rangle _{0,j}}{8\pi m_{\scriptscriptstyle S}^2}\int _{E_\mathrm{{min},i}}^{E_\mathrm{{max},i}} dE \, \frac{dN_{\gamma ,j}}{dE}, \end{aligned}$$

(16)

for an energy bin of width $\varDelta E_i \equiv E_{\text {max},i} - E_{\text {min},i}$, where $\langle \sigma v\rangle _{0,j}\equiv \sigma v_j|_{v\rightarrow 0} \equiv \sigma v_j|_{s\rightarrow 4m_{\scriptscriptstyle S}^2}$ is the partial annihilation cross-section into final state j in the zero-velocity limit, and $dN_{\gamma ,j}/dE$ is the differential photon multiplicity for annihilations into the jth final state.

We use a combination of analytic expressions from Ref. [34] and micrOMEGAs to compute the annihilation and semi-annihilation cross-sections for indirect detection. At tree level, the zero temperature annihilation cross-section for a pair of scalar singlet particles to SM states, $\langle \sigma v\rangle _{0}$, is given by the processes in the top two rows of Fig. 1 for the $\mathbb {Z}_3$ model, and equivalent processes in the $\mathbb {Z}_2$ model with $S=S^*$. The effective cross-sections for annihilation to SM final states in the $\mathbb {Z}_3$ model are a factor of two smaller than in the $\mathbb {Z}_2$ model, accounting for the fact that only particle-antiparticle pairs can annihilate. We compute these by scaling the $\mathbb {Z}_2$ cross-sections down by a factor of two. We obtain the semi-annihilation cross-section directly from micrOMEGAs, with there being no equivalent in the $\mathbb {Z}_2$ model.

With the necessary cross-sections computed we then obtain the predicted spectrum ${dN_{\gamma }}/{dE}$ for each model point by using a Monte-Carlo showering simulation, detailed in Ref. [22]. We then use this to compute a combined likelihood for all the dwarf spheroidals in the Fermi-LAT 6-year Pass 8 dataset [93]. The details of this likelihood are given in Ref. [47].

4.6 Higgs invisible width

When $m_{\scriptscriptstyle S}<m_h/2$, the Higgs may decay to two S bosons (Fig. 1). The resulting S bosons would be invisible at the LHC, so they would be identified as a missing contribution to the total decay width. For a model with a $\mathbb {Z}_3$-charged scalar, the decay width of the Higgs to S bosons is

$$\begin{aligned} \varGamma ^{\mathbb {Z}_3}_{h\rightarrow {\scriptscriptstyle SS^*}} = \frac{\lambda _{h\scriptscriptstyle S}^2 v_0^2}{16\pi m_h}\left( 1 -4 m_{\scriptscriptstyle S}^2/m_h^2\right) ^{1/2}, \end{aligned}$$

(17)

where $v_0$ is the Higgs VEV. In the $\mathbb {Z}_2$ model the final states are identical, so we must include a symmetry factor of 1 / 2 to avoid double counting,

$$\begin{aligned} \varGamma ^{\mathbb {Z}_2}_{h\rightarrow {\scriptscriptstyle SS}} = \frac{1}{2}\varGamma ^{\mathbb {Z}_3}_{h\rightarrow {\scriptscriptstyle SS^*}}\, . \end{aligned}$$

(18)

Equations (17) and (18) show that constraints on the Higgs invisible width exclude large $\lambda _{h\scriptscriptstyle S}$ for small singlet masses.

With SM-like couplings (which the Higgs possesses in the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ models), the upper limit on the invisible branching fraction of the Higgs is 19% at 95% confidence level [95]. We employ the implementation of the full likelihood associated with this result in DecayBit [73].

4.7 Additional likelihoods

We also include simple likelihoods for the nuisance parameters varied in our fits (Table 2), via DarkBit [22] and PrecisionBit [73]. These quantities are well constrained by existing data.

We implement a log-normal likelihood for the local DM density, with a central value of $\bar{\rho _0}=0.4$ GeV/cm$^3$ (e.g. [96]) and an uncertainty of $\sigma _{\rho _0} = 0.15$ GeV cm$^{-3}$,

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

(19)

where $\sigma '_{\rho _0} = \ln (1 + \sigma _{\rho _0}/\rho _0)$. More details can be found in Ref. [22].

Fig. 2

Profile likelihoods for the $\mathbb {Z}_2$ scalar singlet model, with the requirement that $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$ only. Results are shown in the $m_{\scriptscriptstyle S}$–$\lambda _{h\scriptscriptstyle S}$ (top) and $m_{\scriptscriptstyle S}$–$\lambda _{\scriptscriptstyle S}$ (bottom) planes. Left panels show a zoomed-in view of the resonance region; right panels show the full mass range. Contour lines indicate $1\sigma$ and $2\sigma$ confidence regions, and best fit points are indicated with stars. Shading and white contours show the result of including the 2018 XENON1T analysis [55], whereas grey annotations illustrate the impact of using the 2017 analysis [54] instead

We model the speed distribution of DM in the Milky Way as Maxwell–Boltzmann, truncated at the local Galactic escape velocity $v_\mathrm {esc}$. We apply a Gaussian likelihood to the mean of this distribution, characterised by a central value of 240 km s$^{-1}$ and a standard deviation of 8 km s$^{-1}$. This is based on a calculation of the circular rotation speed of the Sun, $v_{\text {rot}}$ [97]. We also constrain the escape velocity using a Gaussian likelihood based on $v_\mathrm {esc} = 550 \pm 35$ km s$^{-1}$, derived from measurements of stellar velocities in the RAVE survey [98].

We apply Gaussian likelihoods to the nuclear parameters as well, based on the estimates $\sigma _s = 43 \pm 8$ MeV [99] and $\sigma _l = 50 \pm 15$ MeV [100]. More detailed discussion of our adopted nuclear and velocity likelihoods can be found in Refs. [22, 101].

For the Higgs mass, the top quark mass and the strong coupling, we use Gaussian likelihoods based on $m_h= 125.09 \pm 0.24$ GeV [71, 102], $m_t = 173.34\pm 0.76$ GeV [71, 103] and $\alpha _{s}(m_{Z}) = 0.1181 \pm 0.0011$ [71].

5 The status of the $\mathbb {Z}_2$ model

5.1 No vacuum constraint

In this section, we present global fits of the $\mathbb {Z}_2$ scalar singlet model, with a full spectrum calculation and RGE running up to the Planck scale. In the most general case, we allow either a metastable or an absolutely stable electroweak vacuum. We furthermore require that the dimensionless couplings remain in the perturbative regime (which we define to be less than $\sqrt{4\pi }$), up to the greater of $\overline{m}_{\scriptscriptstyle S}$ and $m_t$. That is, we demand that $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$. The profile likelihoods for this scenario are presented in Fig. 2 in the $m_{\scriptscriptstyle S}$–$\lambda _{h\scriptscriptstyle S}$ (top) and $m_{\scriptscriptstyle S}$–$\lambda _{\scriptscriptstyle S}$ (bottom) planes. All three of the regions mentioned in Sect. 2 (resonance, neck and high-mass) are clearly visible in the upper panels. Due to the strength of the latest direct detection constraints, and the fact that we rescale the expected signals by the thermal relic density at each point in the parameter space, the high-mass region is split into a TeV-scale mode and an intermediate-mass mode situated just above $m_{\scriptscriptstyle S}=m_h$.⁶

The restriction to $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$ results in a reduction of the volume of the allowed region compared to our results in Ref. [47]. Any model with values of $\lambda _{h\scriptscriptstyle S}$ or $\lambda _{\scriptscriptstyle S}$ greater than $\sqrt{4\pi }$ at the input scale $\overline{m}_{\scriptscriptstyle S}$ violates the perturbativity condition even before RGE running, so our profile likelihoods extend only to $\sqrt{4\pi }$ in $\lambda _{h\scriptscriptstyle S}$ and $\lambda _{\scriptscriptstyle S}$. In contrast, Ref. [47] allowed up to $\lambda _{h\scriptscriptstyle S}= 10$.

At very large $m_{\scriptscriptstyle S}$ and $\lambda _{h\scriptscriptstyle S}$, the Higgs quartic coupling is driven up by large loop corrections in the scalar sector, causing it to become non-perturbative. This excludes a region at $\log _{10}(\lambda _{h\scriptscriptstyle S}) \approx 0.2$, $\log _{10}(m_{\scriptscriptstyle S}/\mathrm {GeV}) \approx 3.6$ that was previously allowed in Ref. [47]. This consequence of the perturbativity requirement at high $m_{\scriptscriptstyle S}$ combines with the relic density constraint, which pushes up at the allowed parameter region from below, to provide a robust upper limit on the singlet mass of $m_s < 4.5\,\mathrm {TeV}$.

Fig. 3

Profile likelihood for the $\mathbb {Z}_2$ scalar singlet model with the requirement $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$. The preferred regions are expressed as a function of $m_{\scriptscriptstyle S}$ and the spin-independent direct detection cross-section for scattering with protons rescaled to the predicted relic abundance $\sigma _p^\text {SI}\cdot f$, and compared to the exclusion bounds from various direct detection experiments. Contour lines indicate $1\sigma$ and $2\sigma$ confidence regions, and best fit points are indicated with stars. Shading and white contours show the result of including the 2018 XENON1T analysis [55], whereas grey annotations illustrate the impact of using the 2017 analysis [54] instead. Coloured solid lines indicate published limits from PandaX [56] and XENON1T [55], and the dashed line is the projected sensitivity of LZ [94]

The profile likelihood of the scalar quartic coupling $\lambda _{\scriptscriptstyle S}$ is reasonably uniform over the prior range. This is unsurprising, given that $\lambda _{\scriptscriptstyle S}$ has little phenomenological impact in this model; indeed, this is why it was not included in Ref. [47]. However, as we will show, it can be important for stabilising the electroweak vacuum and/or influencing the range of scales over which the model can remain perturbative.

In Fig. 3, we show the spin-independent nuclear scattering cross-section $\sigma _p^\text {SI}$ for the $\mathbb {Z}_2$ scalar singlet, rescaled by the fraction $f\equiv \varOmega _{\scriptscriptstyle S} / \varOmega _{\scriptscriptstyle \text {DM}}$ of the relic density explained by each point in parameter space. Because we scale the expected signals of each model by f when computing their likelihoods, this rescaling is necessary when visually comparing predicted cross-sections to published exclusion curves (which assume $f=1$). Compared to our earlier results [47], significant amounts of parameter space are now excluded from the high-mass modes. The new perturbativity constraint removes parameter space at low $\sigma _p^\text {SI}$, whereas the advent of XENON1T cuts into the allowed region from above. LZ [94] will probe a large fraction of the remaining parameter space, even including a substantial part of the resonance region.

In Figs. 2 and 3, we also illustrate the impact of the latest XENON1T data [55] by comparing the 1 and 2$\sigma$ CL regions with the inclusion of the 2018 (white contours) and 2017 constraints (grey contours). We see that the majority of the impact of XENON1T compared to Ref. [47] was provided already in 2017, with a comparatively modest additional constraint imposed by the 2018 data. Indeed, the small excess above background expectation at high recoil energies in the 2018 data leads to a small increase in the size of the 1$\sigma$ preferred region at large $m_{\scriptscriptstyle S}$, where the predictions of the model are consistent with the observed excess.

Fig. 4

Impact of the requirement of vacuum stability on the $\mathbb {Z}_2$ scalar singlet model, expressed in terms of profile likelihoods (left) and posterior probability densities (right). Bullets indicate posterior means and stars indicate best fit points. Shading and white annotations correspond to scans where the singlet is required to absolutely stabilise the electroweak vacuum. For comparison, we also show the preferred regions without this requirement in grey

5.2 Absolutely stable vacuum

Next, we restrict the model further by imposing the additional constraint of absolute vacuum stability (Fig. 4). We find that values of $\lambda _{h\scriptscriptstyle S} > rsim 0.2$ are required to prevent the Higgs quartic coupling becoming negative and thereby stabilise the electroweak vacuum. As a result, the low-mass resonance mode around $m_{\scriptscriptstyle S}\sim m_h/2$ is almost entirely ruled out except for the very top of the neck region, where a few points are found with $\lambda _{h\scriptscriptstyle S} > rsim 0.2$ and a stable vacuum. This essentially leaves just the high mass-modes, centred on approximately $100\,$GeV and 1 TeV, where $\lambda _{h\scriptscriptstyle S}$ is large enough to stabilise the vacuum.

In Fig. 4 we also show the marginalised posterior for the $\mathbb {Z}_2$ model. As in Ref. [47], we see that even without the requirement of absolute vacuum stability, there is a clear preference for the high-mass region over the resonance region, due to the need to fine-tune nuisance parameters in order to fit all existing data at any given point in the resonance region. With the inclusion of vacuum stability, the same effect can be seen to disfavour the medium-mass mode, where $m_{\scriptscriptstyle S}$ is $\mathcal {O}(100)$ GeV.

Both the profile likelihoods and the marginalised posterior of Figs. 2 and 4 show a small diagonal strip where valid solutions are difficult to come by at large $m_{\scriptscriptstyle S}$ and $\lambda _{h\scriptscriptstyle S}$, just below (and running parallel to) the border of the allowed region where the Higgs quartic coupling becomes non-perturbative. In this region, the pole mass calculation for the Higgs runs into numerical instabilities, and fails to converge. This is a numerical artefact; large $\lambda _{h\scriptscriptstyle S}$-dependent radiative corrections cause the Higgs pole mass calculation at $\overline{m}_{\scriptscriptstyle S}$ to fail in the UV (singlet) theory, but technically this particular iteration could be avoided, seeing as the Higgs pole mass that we actually adopt comes instead from the SM EFT. The true results in this region would therefore smoothly interpolate those from the surrounding region. This effect confirms that predictions are becoming less stable, due to large one-loop corrections, as we approach the perturbativity limit, and that indeed we should not adopt any larger perturbative cutoff on the dimensionless couplings than $\sqrt{4\pi }$. This problem can also be partially compensated for by varying other nuisance parameters (in particular, the top mass) within their allowed ranges, as can be seen by the fact that this effect has a much larger impact on the posterior than the profile likelihood.

Fig. 5

Scale of perturbativity violation with respect to $m_{\scriptscriptstyle S}$ and $\lambda _{\scriptscriptstyle S}$ for the $\mathbb {Z}_2$ scalar singlet model with the requirement that $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$ only (left) and with the additional requirement of absolute vacuum stability (right). The $1\sigma$ and $2\sigma$ confidence regions are delineated by white contours, and the best-fit by a white star

Let us now take a closer look at just how the $\mathbb {Z}_2$ scalar singlet model can satisfy the vacuum stability constraint. As discussed in Sect. 4.2, we apply the vacuum stability condition by excluding points where we can show perturbatively that the electroweak vacuum is metastable, because $\lambda _{h}$ becomes negative before the Planck scale $M_{\text {Pl}}=1.22\times 10^{19}$ GeV. However, this means that in Fig. 4, we do not distinguish between two quite different cases:

• (i) at high scales all couplings remain perturbative and $\lambda _{h}\ge 0$,

• (ii) some couplings simply run to non-perturbative values before $M_{\text {Pl}}$.

In the case of (i), we have explicitly shown that the scalar singlet model can help to stabilise the electroweak vacuum. In (ii), the stability of the electroweak vacuum may be restored by non-perturbative effects, but we are unable to determine whether this is the case or not from our perturbative calculations. It is therefore important to discriminate between these two cases.

In Fig. 5, we plot the scale at which perturbativity is violated in the $m_{\scriptscriptstyle S}$–$\lambda _{\scriptscriptstyle S}$ parameter plane, choosing the scale by profiling the likelihood over the other parameters (i.e. plotting $\varLambda _P$ for the best-fit points found in the scan at each combination of $m_{\scriptscriptstyle S}$ and $\lambda _{h\scriptscriptstyle S}$). We plot the value of $\varLambda _P$ only within the $2\sigma$ contours, as determined by the profile likelihood. Note that there can exist parameters with a larger value of $\varLambda _P$ that only have a slightly worse $\mathcal {L}$, and are still within $2\sigma$ of the best-fit point. As we run the couplings to a maximum scale of $1\times 10^{20}\,$GeV (well above the Planck scale, where quantum gravitational effects become important), points with $\varLambda _P$ equal to this value should be interpreted as valid to at least $1\times 10^{20}\,$GeV.

Fig. 6

Left: Scale of perturbativity violation for the $\mathbb {Z}_2$ scalar singlet model with the requirement of a stable electroweak vacuum. Right: Profile likelihood when furthermore imposing the requirement $\varLambda _P>10^{15}\,$GeV. The $1\sigma$ and $2\sigma$ confidence regions are delineated by white contours, and the best-fit by a white star. Grey contours on the right panel correspond to the $1\sigma$ and $2\sigma$ confidence regions of the left panel

A number of observations can be made from Fig. 5. First of all, we note that for $\lambda _{\scriptscriptstyle S} > rsim 0.7$ the scale of perturbativity is always low, as the singlet quartic coupling quickly runs to non-perturbative values. The dependence on $m_{\scriptscriptstyle S}$ is more complicated and can be best understood by comparing to Fig. 2. In the low-mass resonance mode ($m_{\scriptscriptstyle S}\sim m_h/2$), $\lambda _{h\scriptscriptstyle S}$ is typically very small and the scale of perturbativity violation can be very high as long as $\lambda _{\scriptscriptstyle S}$ is sufficiently small. The mode at $m_{\scriptscriptstyle S}\sim 100\,$GeV, on the other hand, requires $\lambda _{h\scriptscriptstyle S}> 1$, which renders the spectrum invalid at scales well below $10^{10}\,$GeV irrespective of $\lambda _{\scriptscriptstyle S}$. In the high-mass region ($m_{\scriptscriptstyle S}\sim 1\,\mathrm {TeV}$) it is possible to find points with a scale of perturbativity near or beyond the Planck scale, in particular towards smaller masses (corresponding to smaller $\lambda _{h\scriptscriptstyle S}$).

In Fig. 6 (left), we show $\varLambda _P$ as function of $\lambda _{h\scriptscriptstyle S}$ and $m_{\scriptscriptstyle S}$ in the high-mass region, imposing absolute vacuum stability. There is a rough correlation between $\lambda _{h\scriptscriptstyle S}$ and $\varLambda _P$, which is only broken for the small $\lambda _{h\scriptscriptstyle S}$ tip of the high-mass mode, where $\varLambda _P$ decreases rapidly. This is because such small values of $\lambda _{h\scriptscriptstyle S}$ are insufficient to stabilise the electroweak vacuum, so our requirement that $\lambda _{h}$ not run negative only finds solutions where $\lambda _{\scriptscriptstyle S}$ contributes to the running of $\lambda _{h}$. However, since the impact of $\lambda _{\scriptscriptstyle S}$ on the running of $\lambda _{h}$ is indirect and only mild, large values of $\lambda _{\scriptscriptstyle S}$ are required, rendering the model non-perturbative below the scale of vacuum instability (thus rendering it “stable” according to our definition). This can also be seen as the cause for the difference in the profile likelihood and the posterior in the tip of this region in Fig. 4, reflecting the fact that $\lambda _{\scriptscriptstyle S}$ must be tuned in order to find permitted models in this area.

The competing interests of vacuum stability and perturbativity become more problematic when we ask what values of $\varLambda _P$ are acceptable. The metastability of the electroweak vacuum in the SM is the result of the Higgs quartic coupling becoming negative near the grand unified theory (GUT) scale, at $\sim \,10^{15}\,$GeV. If we are concerned about vacuum stability, then we should generally also demand that our theory is perturbative to at least this scale. The electroweak vacuum is stable in the $\mathbb {Z}_2$ theory in some parts of the otherwise allowed parameter space, but in others the model simply becomes non-perturbative at scales well below the GUT scale. It therefore makes sense to impose another selection requirement on our samples, in order to identify only those points that remain perturbative to high scales.

In the right panel of Fig. 6, we show the profile likelihood after excluding all models where $\varLambda _P<10^{15}\,$GeV (in addition to requiring a stable electroweak vacuum). Due to the competing influences of direct detection, the relic density, vacuum stability and perturbativity, the allowed parameter space of the model is reduced considerably. Nevertheless, even when limiting the parameter space to points with $\varLambda _P>10^{15}\,$GeV, the model is certainly not ruled out. In the next section, we will explicitly identify parameter points for which couplings remain perturbative up to the typical instability scales, and the electroweak vacuum is stabilised, showing that these points can still give a good fit to the data – and can even explain the entirety of DM.

In Fig. 7, we show the impacts of demanding absolute vacuum stability and perturbativity to high scales on signals at direct detection experiments. If the $\mathbb {Z}_2$ model is to stabilise the electroweak vacuum, it must lie in a narrow region with effective nuclear scattering cross-section $5\times 10^{-46}< \sigma _p^\text {SI}\cdot f < 10^{-45}$ cm$^2$ and mass 600 GeV < $m_{\scriptscriptstyle S}< 2$ TeV. Here we again show both the result with the 2018 (shading and white contours) and 2017 XENON1T likelihoods (grey contours). Even with vacuum stability imposed, the new XENON1T data remains consistent with the allowed region. Future multi-ton experiments such as LZ [94] will definitively detect or exclude this scenario.

Fig. 7

As in Fig. 3, but with the added requirements of vacuum stability and perturbativity to large scales, $\varLambda _P>10^{15}\,$GeV

Table 5 Details of the best-fit points for the $\mathbb {Z}_2$ scalar singlet model when different physical restrictions are imposed on the model. Points that have an absolutely stable electroweak vacuum are indicated by a tick in the first column. Points with a singlet relic density within 1$\sigma$ of the Planck observed value ($\varOmega _{\scriptscriptstyle S} h^2\sim \varOmega _{\scriptscriptstyle \text {DM}} h^2$) are indicated with a tick in the third column. We omit the values of the nuisance parameters, as they are not significantly different to the central values of their respective likelihood functions

Table 6 Individual contributions to the $\varDelta$ log-likelihood for the various best-fit points (see Tables 5 and 7) compared to an ‘ideal’ case. We take this to be the background-only likelihood for exclusions, and the central observed value for detections. Note that because each likelihood is dimensionful, the absolute values are less meaningful than the offset with respect to another point (see section 8.3 of Ref. [15] for more details on the normalisation). The best-fit points are labelled as follows: A represents a fit with the only constraint that $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$, B is a fit with the additional constraint of absolute vacuum stability, C includes the constraint of $\varLambda _P>10^{15}$  GeV and D also includes the requirement that $\varOmega _Sh^2$ be within 1$\sigma$ of the observed relic density

5.3 Best-fit point

The best-fit point for our global fit of the $\mathbb {Z}_2$ model is located at $\lambda _{\scriptscriptstyle S}= {6.24\times 10^{-3}}$, $\lambda _{h\scriptscriptstyle S}= {2.32\times 10^{-4}}$ and $m_{\scriptscriptstyle S}= {62.48}\,$GeV. This point is located in the low-mass resonance region, the electroweak vacuum is metastable with a lifetime of $\sim$ ${1.1\times 10^{99}}$ years, the minimum of $\lambda _{h}$ occurs at $\sim 3\times 10^{13}\,$GeV, and the model is perturbative up to at least $10^{20}\,$GeV. Details of this point can be found in Table 5. The mass at this point is within $0.03\,$GeV of the best-fit found in Ref. [47], and the portal coupling is approximately a factor of three smaller. Given that the profile likelihood is quite flat with respect to $\lambda _{h\scriptscriptstyle S}$ around the best-fit, the difference in $\lambda _{h\scriptscriptstyle S}$ is not significant.

Indeed, we expect this best-fit point to be similar to what we found in Ref. [47], as the constraint $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$ and the variation of $\lambda _{\scriptscriptstyle S}$ do not have a significant impact on the phenomenology at small couplings. Nevertheless, we find $\varDelta \ln \mathcal {L} = 0.317$ relative to the ideal likelihood (where each individual likelihood takes its maximum value), compared to $\varDelta \ln \mathcal {L} = 0.107$ in Ref. [47]. This difference is due to the contribution from the new 2018 XENON1T likelihood, which exhibits a slight preference for a non-zero DM signal and hence slightly disfavours the low mass region (see Table 6).

When the constraint of absolute vacuum stability is imposed, the location of the best fit necessarily moves away from the resonance region, where $\lambda _{h\scriptscriptstyle S}$ is too small to stabilise the vacuum. In this case we find a best-fit point at $\lambda _{\scriptscriptstyle S}={2.28\times 10^{-3}}$, $\lambda _{h\scriptscriptstyle S}={2.03}$ and $m_{\scriptscriptstyle S}={3.97}\,$TeV. In this case we find $\varDelta \ln \mathcal {L} = 0.503$, which corresponds to a slight penalty over the metastable case. Note in particular that this second point gives a better fit to the data from XENON1T, but the combined likelihood from all direct detection experiments is worse due to the contribution from LUX and PandaX. Although the vacuum is classified as stable at this point, it is an example of a point where the couplings are so large that they cannot be run all the way to the typical scale of vacuum instability, leading to $\varLambda _P\sim 26\,$TeV. This reduces the theoretical appeal of this point.

Fig. 8

Profile likelihoods for the $\mathbb {Z}_3$ scalar singlet model, with the requirement that $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$ only. Results are shown in the $\lambda _{h\scriptscriptstyle S}$–$m_{\scriptscriptstyle S}$ (left) and in the $\mu _3$–$m_{\scriptscriptstyle S}$ (right) planes. Contour lines indicate $1\sigma$ and $2\sigma$ confidence regions, and best fit points are indicated with stars. Shading and white contours show the result of including the 2018 XENON1T analysis [55], whereas grey annotations illustrate the impact of using the 2017 analysis [54] instead

By excluding all samples with $\varLambda _P<10^{15}\,$GeV, we can find points that have a stable vacuum and are more theoretically interesting (see Fig. 6). We find a best-fit point that is absolutely stable and has $\varLambda _P = 1.0\times 10^{15}\,$GeV, at $\lambda _{\scriptscriptstyle S}={6.59\times 10^{-4}}$, $\lambda _{h\scriptscriptstyle S}= {7.36\times 10^{-1}}$, $m_{\scriptscriptstyle S}= {1.93}\,$TeV. This point has $\varDelta \ln \mathcal {L} = 1.121$, with the largest contributions coming from direct detection likelihoods. This corresponds to a likelihood ratio $\varLambda = 0.448$ relative to the overall best-fit point, which places it inside the generally-preferred $1\sigma$ parameter region.

Although the requirement of an absolutely stable vacuum and perturbative couplings up to at least the GUT scale leads to some mild tension with direct detection limits, it is intriguing to observe that this tension has not grown with the latest XENON1T results, even though the expected sensitivity of XENON1T would have been sufficient to comprehensively test these solutions. The reason is that in precisely this parameter region, the model accommodates the slight preference for a non-zero DM signal in XENON1T. It will therefore be extremely interesting to include the results from the next generation of direct detection experiments in a similar analysis.

Finally, we consider points with relic densities within $1\sigma$ of the Planck measured value, with stable vacua, and that remain perturbative to at least $10^{15}\,$GeV. The best-fitting of these points is located at $\lambda _{\scriptscriptstyle S}= {2.59\times 10^{-3}}$, $\lambda _{h\scriptscriptstyle S}= {6.80\times 10^{-1}}$ and $m_{\scriptscriptstyle S}= {1.94}\,$TeV, and has $\varDelta \ln \mathcal {L} = 1.455$, still within 1$\sigma$ of the global best-fit point. The four best-fit points, the corresponding relic densities and the scale of perturbativity violation are presented in Table 5. The individual likelihood contributions for each point are given in Table 6.

By interpreting $\varDelta \ln \mathcal {L}$ as half the “likelihood $\chi ^2$” of Baker and Cousins [104] and assuming either one or two degrees of freedom, we can obtain an approximate p-value for each of our best-fit points. For the best fit with metastability allowed, we find $p\approx 0.4\text {--}0.7$. With vacuum stability required, p drops to 0.3–0.6. For the case where the couplings are perturbative up to $10^{15}\,$GeV and the electroweak vacuum is absolutely stable, we find $p \approx 0.15\text {--}0.3$. This decreases to $p \approx 0.1\text {--}0.25$ when also requiring that the S relic density is within $1\sigma$ of the Planck value. Each of these p-values is acceptable, although requiring the UV properties of perturbativity and vacuum stability does have a notable impact.

6 The status of the $\mathbb {Z}_3$ model

6.1 No vacuum constraint

We now turn to the $\mathbb {Z}_3$ scalar singlet model. The main difference compared to the $\mathbb {Z}_2$ model is the presence of an additional parameter $\mu _3$, which has a significant impact on phenomenology because it can lead to semi-annihilations. Fig. 8 presents the profile likelihoods in the $m_{\scriptscriptstyle S}$–$\lambda _{h\scriptscriptstyle S}$ parameter plane (left) and in the $m_{\scriptscriptstyle S}$–$\mu _3$ parameter plane (right), based on scans over the full range of $m_{\scriptscriptstyle S}$. Note that the allowed region for $\mu _3$ is constrained by the vacuum stability condition given in Eq. (11), particularly at small singlet masses.

Fig. 9

The semi-annihilation fraction $\alpha$ with respect to $\lambda _{h\scriptscriptstyle S}$ and $m_{\scriptscriptstyle S}$ (left) and with respect to $\mu _3$ and $m_{\scriptscriptstyle S}$ (right) for the $\mathbb {Z}_3$ scalar singlet model without the requirement of an absolutely stable electroweak vacuum, but imposing $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$. The $1\sigma$ and $2\sigma$ confidence regions are delineated by white contours, and the best-fit by a white star

Fig. 10

Profile likelihood for the $\mathbb {Z}_3$ scalar singlet model with the requirement $\varLambda _P>\max (\overline{m}_{\scriptscriptstyle S},m_t)$. Regions are shown as a function of $m_{\scriptscriptstyle S}$ and the spin-independent direct detection cross-section for scattering with protons, rescaled by the predicted relic abundance $\sigma _p^\text {SI}\cdot f$, and compared to the exclusion bounds from various direct detection experiments. Contour lines indicate $1\sigma$ and $2\sigma$ confidence regions, and best fit points are indicated with stars. Shading and white contours show the result of including the 2018 XENON1T analysis [55], whereas grey annotations illustrate the impact of using the 2017 analysis [54] instead. Other lines indicate limits from PandaX [56] and XENON1T [55], and the projected sensitivity of LZ [94]

We find that the allowed resonance region in the $\mathbb {Z}_3$ model is practically identical to the corresponding region in the $\mathbb {Z}_2$ model. We therefore do not include a version of Fig. 8 zoomed in to low masses. At larger masses, however, there are notable differences between the $\mathbb {Z}_2$ and $\mathbb {Z}_3$ models. The allowed parameter region with $m_{\scriptscriptstyle S}\sim 200\,$GeV is substantially larger in the $\mathbb {Z}_3$ model, and extends to much smaller values of $\lambda _{h\scriptscriptstyle S}$. This difference in shape can be understood by considering the fraction of semi-annihilation. In Fig. 9 we plot the semi-annihilation fraction $\alpha$, defined in Eq. (14), within the $2\sigma$ confidence regions. As expected, the extended allowed parameter region in the intermediate mass range corresponds to $\alpha \approx 1$, meaning that the semi-annihilation channel dominates. As a result, the same relic abundance can be achieved with smaller values of the portal coupling $\lambda _{h\scriptscriptstyle S}$. In other words, the bound $\varOmega _{\scriptscriptstyle S}h^2\le \varOmega _{\scriptscriptstyle \text {DM}}h^2$ can be evaded at much lower values of $\lambda _{h\scriptscriptstyle S}$, by invoking a large contribution from semi-annihilation. Similarly, at large values of $\lambda _{h\scriptscriptstyle S}$ the contribution from semi-annihilation brings the relic density lower than in the $\mathbb {Z}_2$ model, and direct detection constraints are more easily avoided (given that we consistently rescale signals for the local density of singlet particles).

Fig. 11

Impact of the requirement of vacuum stability on the $\mathbb {Z}_3$ scalar singlet model, expressed in terms of profile likelihoods (top) and posterior probability densities with flat prior on $\mu _3$ (centre) and with logarithmic prior on $\mu _3$ (bottom). Bullets indicate posterior means and stars indicate best fit points. Shading and white annotations correspond to scans where the singlet is required to absolutely stabilise the electroweak vacuum. For comparison, we also show the preferred regions without this requirement in grey

For larger masses, semi-annihilations become less efficient, as the semi-annihilation fraction is proportional to $\mu _3\lambda _{h\scriptscriptstyle S}^2/m_{\scriptscriptstyle S}^6$ at leading order [8]. As a result, the shape of the allowed parameter region is similar to the $\mathbb {Z}_2$ model for $m_{\scriptscriptstyle S} > rsim 1 \, \mathrm {TeV}$. The likelihood of this region is however much smaller, and is in fact outside the global $1\sigma$ confidence region. The reason is that the $\mathbb {Z}_3$ model requires a complex scalar, whereas we have considered a real scalar for the $\mathbb {Z}_2$ model. The coupling $\lambda _{h\scriptscriptstyle S}$ must therefore be a factor of two larger in the $\mathbb {Z}_3$ model to achieve the necessary effective annihilation cross-section required to avoid DM overproduction, increasing the degree of tension between the relic density constraint and bounds from direct detection experiments.

Fig. 12

Left: Scale of perturbativity violation for the $\mathbb {Z}_3$ scalar singlet model with the requirement of a stable electroweak vacuum. Right Profile likelihood when also imposing the requirement $\varLambda _P>10^{15}\,$GeV. The $1\sigma$ and $2\sigma$ confidence regions are delineated by white contours, and the best-fit by a white star. Grey contours on the right panel correspond to the $1\sigma$ and $2\sigma$ confidence regions of the left panel

Figure 10 illustrates the impact of current and future direct detection experiments on the parameter space of the $\mathbb {Z}_3$ model. Here we have again rescaled $\sigma _p^\text {SI}$ by the relic density fraction $f = \varOmega _s/\varOmega _\text {DM}$. Similar to the case of the $\mathbb {Z}_2$ model, the resonance region extends three orders of magnitude below current direct detection limits, and even a next-generation experiment such as LZ will not be able to probe the full parameter space at $m_{\scriptscriptstyle S}\simeq m_h/2$. However, the allowed parameter region at $m_{\scriptscriptstyle S}\sim 200\,$GeV, even though it is substantially larger than its counterpart in the $\mathbb {Z}_2$ model, will eventually be fully explored by direct detection searches.

6.2 Absolutely stable vacuum

As with the $\mathbb {Z}_2$ model, we now investigate the preferred parameter regions more closely by imposing additional physical requirements. In the top panel of Fig. 11, we show how the profile likelihoods change when requiring absolute vacuum stability. The corresponding marginalised posteriors are shown in the central panel of Fig. 11 for a scan with flat prior on $\mu _3$, and in the bottom panel for a scan with logarithmic prior on $\mu _3$. Although the likelihood is maximised for $m_{\scriptscriptstyle S}\sim 100\,\mathrm {GeV}$, the majority of the posterior mass is found at higher masses, $m_{\scriptscriptstyle S}> 1\,\mathrm {TeV}$. The reason is that at larger masses, the relic abundance becomes independent of $\mu _3$ and therefore benefits strongly from marginalisation (rather than profiling) over $\mu _3$. Comparing the middle and lower panels, this conclusion is independent of the choice of prior on $\mu _3$. The choice of prior on $\mu _3$ has a relatively small impact overall, essentially just translating into a stronger preference for large $m_{\scriptscriptstyle S}$ when taken flat rather than logarithmic, due to the restriction to lower values of $\mu _3$ at lower $m_{\scriptscriptstyle S}$ coming from Eq. 11.

Figure 11 demonstrates that most of the parameter space opened up by semi-annihilations in the intermediate mass range remains viable when imposing absolute vacuum stability. This observation raises the question whether the stabilisation is due to the influence of $\mu _3$ on the running of $\lambda _{h}$ or whether the new parameter simply leads to a breakdown of perturbativity. We therefore show the scale of perturbativity violation in the left panel of Fig. 12. Indeed, we find that $\varLambda _P$ is extremely low (less than $10^{10}\,$GeV) throughout the $2\sigma$ preferred region of the new parameter space, such the points with a ‘stable’ vacuum are not actually as theoretically appealing as might have naively been expected on the basis of Fig. 11. This is also the reason for the persistence of part of the resonance reason in Fig. 11 after absolute vacuum stability is required, unlike in the corresponding $\mathbb {Z}_2$ plot (Fig. 4). $\mathbb {Z}_3$ models remaining in this region after absolute vacuum stability is required are just those with the highest values of $\lambda _{\scriptscriptstyle S}$, allowing them to combine with non-zero values of $\mu _3$ to send $\lambda _{\scriptscriptstyle S}$ non-perturbative at relatively low scales, and thereby avoid having to actually stabilise the electroweak vacuum.

The reason that large couplings are required in the intermediate mass range is related to the need for semi-annihilations in this part of parameter space. A large semi-annihilation fraction requires that the coupling $\mu _3$ is larger than about $300\,$GeV (see Fig. 9). This in turn forces $\lambda _{\scriptscriptstyle S}$ to be large in order to satisfy Eq. (11), which leads to the couplings becoming non-perturbative at a low scale.⁷ Thus, the new regions of parameter space opened up by semi-annihilations in the $\mathbb {Z}_3$ model are of limited theoretical appeal from the point of view of stabilising the electroweak vacuum.

Fig. 13

As in Fig. 10, but with the added requirements of vacuum stability and perturbativity to large scales, $\varLambda _P>10^{15}\,$GeV

Table 7 Details of the best-fit points for the $\mathbb {Z}_3$ scalar singlet model when different physical restrictions are imposed on the model. Points that have an absolutely stable electroweak vacuum are indicated by a tick in the first column. Points with a singlet relic density within 1$\sigma$ of the Planck observed value ($\varOmega _{\scriptscriptstyle S} h^2\sim \varOmega _{\scriptscriptstyle \text {DM}} h^2$) are indicated with a tick in the third column. We omit the values of the nuisance parameters, as they are not significantly different to the central values of their respective likelihood functions

In the right panel of Fig. 12, we present the profile likelihood after also imposing the requirement $\varLambda >10^{15}\,$GeV. We find that this not only removes the resonance region,⁸ but also the intermediate mass range. In fact, only a tiny region around $m_{\scriptscriptstyle S}\sim 1\,\mathrm {TeV}$ remains. Judging from Fig. 8, we expect this remaining parameter space to have a much lower likelihood than the resonance region. In Fig. 13 we show the nuclear scattering cross-section (rescaled as usual by the fraction f of the DM relic density constituted by singlet scalars) as a function of $m_{\scriptscriptstyle S}$. This plot qualitatively confirms our expectation that this region should lie outside the preferred region in the global scan, as the 90% C.L. upper bound from XENON1T [55] already excludes the $2\sigma$ confidence regions. In the following, we will make this point more explicit by studying the best-fit point in this region, and showing that it is in considerable tension with data.

6.3 Best-fit point

The best-fit point in the $\mathbb {Z}_3$ model with metastability allowed is at $\lambda _{\scriptscriptstyle S}={4.83\times 10^{-1}}$, $\lambda _{h\scriptscriptstyle S}={3.21\times 10^{-4}}$, $\mu _3 = 11.8\,$GeV and $m_{\scriptscriptstyle S}= {62.48}\,$GeV. This point has a lifetime of $\sim {1.4\times 10^{99}}$ years, and a minimum in its Higgs quartic coupling at $\sim 3\times 10^{13}\,$GeV. For this point we find $\varDelta \ln \mathcal {L} = 0.318$, essentially the same as for the equivalent best-fit point in the $\mathbb {Z}_2$ model. This result is expected, as the semi-annihilation fraction at this point is $\alpha = 0$.

With the additional constraint of absolute vacuum stability, the best-fit is located at $\lambda _{\scriptscriptstyle S}={3.52}$, $\lambda _{h\scriptscriptstyle S}={3.50}$, $\mu _3 = {537}\,$GeV and $m_{\scriptscriptstyle S}={144}\,$GeV. In this case the semi-annihilation fraction is $\alpha =0.72$ and we find $\varDelta \ln \mathcal {L} = 0.473$. Compared to the equivalent point in the $\mathbb {Z}_2$ model, this represents a small improvement due to the contribution from semi-annihilations. However, $\varLambda _P$ at this point is only ${168}\,$GeV, due to the large value of $\lambda _{\scriptscriptstyle S}$, making it less than appealing in a theoretical sense.

Demanding that $\varLambda _P\ge 10^{15}\,$GeV, we find a best-fit point with an absolutely stable vacuum, $\varLambda _P = 1.29\times 10^{15}\,$GeV and $\alpha =0.004$. This point is located at $\lambda _{\scriptscriptstyle S}= {3.08\times 10^{-2}}$, $\lambda _{h\scriptscriptstyle S}= {6.60\times 10^{-1}}$, $\mu _3 = {342}\,$GeV and $m_{\scriptscriptstyle S}= {1.31}\,$TeV. This point has $\varDelta \ln \mathcal {L} = 3.975$, with the dominant contributions coming from the most recent direct detection experiments. This corresponds to a likelihood ratio $\varLambda =0.026$, which places this point more than $2\sigma$ away from the overall best-fit point. In other words, we find considerable tension in the $\mathbb {Z}_3$ model between direct detection limits and the requirement for the model to be absolutely stable and perturbative to at least the GUT scale.

Finally, we consider a point that is perturbative to at least $10^{15}\,$GeV, has a stable electroweak vacuum and has a singlet relic density within $1\sigma$ of the Planck measured value. The best-fit point under these requirements is located at $\lambda _{\scriptscriptstyle S}={1.19\times 10^{-1}}$, $\lambda _{h\scriptscriptstyle S}= {5.92\times 10^{-1}}$, $\mu _3 = 540\,$GeV and $m_{\scriptscriptstyle S}= {1.21}\,$TeV. This point has $\varOmega _Sh^2=0.1137$ and $\varLambda _P = 1.89\times 10^{15}\,$GeV as well as $\varDelta \ln \mathcal {L} = 4.443$, making it even more strongly disfavoured than the corresponding best-fit model with subdominant singlet DM.

We present the four best-fit points, their relic densities and the scales at which their couplings become non-perturbative in Table 7.

As with the $\mathbb {Z}_2$ model, we can obtain approximate p-value ranges by assuming either one or two degrees of freedom. For the best-fit point with metastability allowed, we find $p\approx 0.4\text {--}0.7$, while the best fit with vacuum stability has $p\approx 0.3\text {--}0.6$. Both of these ranges are very similar to those for equivalent constraints on the $\mathbb {Z}_2$ model, despite semi-annihilations opening up a large region of parameter space. For the model with $\varLambda _P > 10^{15}\,$GeV, we find $p\approx 0.005\text {--}0.02$, reducing to $p\approx 0.003\text {--}0.01$ when also imposing the relic density requirement. This illustrates once again that the $\mathbb {Z}_3$ model, which requires a complex scalar, is disfavoured by data as a joint mechanism to stabilise the electroweak vacuum and to provide a DM candidate.

7 Conclusions

In this work we have investigated two realisations of scalar singlet DM: a real scalar stabilised by a $\mathbb {Z}_2$ symmetry and a complex scalar stabilised by a $\mathbb {Z}_3$ symmetry. In addition to potentially accounting for the observed DM relic abundance via the freeze-out mechanism, these models have the attractive feature that in large regions of parameter space they remain valid up to very large scales. This makes it possible to study the RGE evolution of the various parameters and determine the impact of the scalar singlet on the running of the Higgs quartic self-coupling. Indeed, we find that this additional contribution may stabilise the electroweak vacuum, thus resolving an apparent deficiency of the SM.

Nevertheless, models of scalar singlets face a large number of experimental and theoretical constraints. The most important experiments are those aimed at direct detection of DM, most notably the very recent results from XENON1T [55]. In spite of observing a small upward fluctuation, XENON1T places strong constraints on the DM-nucleon scattering cross-section. The most important theoretical requirement that we have considered is for couplings to remain perturbative up to the scales where the stability of the electroweak vacuum may become an issue.

By performing a global fit to all available data, we have shown that it is still possible to explain DM and stabilise the electroweak vacuum through the addition of a scalar singlet field charged under a $\mathbb {Z}_2$ symmetry, while at the same time satisfying all experimental constraints. Although in much of the allowed parameter space we find the scale at which couplings become non-perturbative to be quite low, there is an allowed parameter region with scalar masses of about $2\,$TeV, where the theory remains perturbative up to at least $10^{15}\,$GeV. Moreover, it is possible in this parameter region for scalar singlets to constitute all of the DM. This parameter region is in slight tension with direct detection, but is consistent with the small preference for a non-zero signal contribution at high DM mass in recent XENON1T results. The next generation of experiments will therefore fully explore the viability of this scenario.

The alternative possibility of a complex scalar singlet with a $\mathbb {Z}_3$ symmetry opens up large regions of the parameter space, because the semi-annihilation channel allows the same relic density to be achieved for much smaller singlet-Higgs couplings than in equivalent parts of the $\mathbb {Z}_2$ parameter space. However, the presence of a large trilinear coupling drives the couplings non-perturbative at a relatively low scale, making it impossible to calculate the running of Higgs self-couplings to high scales. When requiring a stable electroweak vacuum as well as perturbative couplings up to at least $10^{15}\,$GeV, the semi-annihilation channel ceases to be relevant and the remaining parameter space resembles the one of the $\mathbb {Z}_2$ model. However, the relic density constraint is more severe for a complex scalar than for a real scalar, so that the $\mathbb {Z}_3$ model is in fact more tightly constrained. Indeed, the parameter region with $m_{\scriptscriptstyle S}\sim 1\,$TeV is disfavoured more than 95% confidence, irrespective of whether or not the scalar singlets constitute all of DM.

Scalar singlets have frequently been advocated as one of the simplest realisations of the WIMP idea. The non-observation of a DM signal in any type of experiment designed to search for WIMPs therefore clearly increases pressure on these models. While the low-mass (resonance) region remains challenging to probe experimentally, we have focused on the more interesting high-mass region, which may help to address the issue of a metastable electroweak vacuum. While this solution is now essentially ruled out for the case of a complex scalar singlet (stabilised e.g. by a $\mathbb {Z}_3$ symmetry), it remains an interesting possibility for real scalars (with a $\mathbb {Z}_2$ stabilising symmetry). The next generation of direct detection experiments will be able to reach a definite verdict on these models, including the exciting possibility that they may confirm the slight excess seen in XENON1T.

Appendix A: Supplementary material

All samples, input files and best-fit points from this paper are freely available for download from Zenodo [57].