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

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 to 2 TeV and a spin-independent nuclear scattering cross-section around 10$^{-45}$ cm$^2$.


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 a j.mckay14@imperial.ac.uk b p.scott@imperial.ac.uk 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 Pl ∼ 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 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].
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 [12], which enables the user to incorpo-rate existing software via a backend system. Specifically, we use FlexibleSUSY 2.0.1 [13,14] and SARAH 4.12.2 [15][16][17][18] for the renormalisation group evolution needed to study vacuum stability, DDCalc 2.0.0 [19] for DM direct detection, gamLike 1.0.0 [19] and DarkSUSY 5.1.3 [20,21] for DM indirect detection with gamma-rays, mi-crOMEGAs 3.6.9.2 [22] for the relic density calculation, and Diver 1.0.4 and T-Walk 1.0.1 [23] 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 Z 2 symmetry stabilises the singlet, and makes it a viable DM candidate . Perturbativity, vacuum stability and the relic density of scalar singlets with a Z 3 symmetry have also been investigated [6]. 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 [6], and impact indirect detection by modifying the injection spectra of light particles [46]. However, vacuum stability considerations limit the magnitude of the responsible coupling, leading to an interesting interplay of different constraints.
In this paper, we present an updated global analysis of the Z 2 model and carry out the first global fit of the Z 3 model. In particular, we improve on our earlier analysis of the Z 2 model [44] by treating the singlet self-coupling λ S as a free parameter, including full RGE running, considering vacuum stability and perturbativity, and incorporating the latest results from XENON1T [47,48] and PandaX [49]. 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 Sec. 2, of our input parameters and scanning procedure in Sec. 3, and of our observable calculations and likelihood functions in Sec. 4. The results for the Z 2 and Z 3 models appear in Secs. 5 and 6, respectively. We summarise our findings in Sec. 7.
GAMBIT software can be downloaded from gambit.hepforge.org, and all samples, input files and best-fit points from this paper are available from Zenodo [50].

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 Z 2 transformation S → −S. The most general renormalisable scalar potential permitted by the Z 2 , Lorentz and gauge symmetries is then [51] The terms proportional to µ 2 H and µ 2 S are the Higgs and singlet bare masses, the terms proportional to λ h and λ 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 µ 2 S , λ hS and λ S . After electroweak symmetry breaking we can replace H → 0, h being the SM Higgs field and v 0 = 246 GeV the VEV of the electroweak vacuum. Then, the portal term proportional to λ hS induces couplings of h to the scalar singlet S via the terms h 2 S 2 and v 0 hS 2 . Moreover, after symmetry breaking the M S singlet mass is given by where v 0 is the M S Higgs VEV. The singlet pole mass, m S , can be obtained from this using where Σ S represents loop corrections that shift the M S mass to the pole. The only renormalisable interaction of S with the SM is through the "Higgs portal" S 2 H 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 [25][26][27], direct detection and h → SS decays [28]. Notice that for scalar masses less than a few TeV, the couplings λ S and λ hS 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 [9,52,53] and baryogenesis [54][55][56].
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 S and λ hS . Specifically, the parameters that have been identified to be compatible with current experimental data exist in a number of regions [31,43,44]: spite of very small couplings (λ hS 10 −2 ) the singlet can nevertheless account for the entire observed relic abundancce of DM. 2. The resonant "neck" region at m 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 λ hS of order one.
Eventually, direct detection is expected to probe much of this remaining parameter space, leaving only large values of λ hS at which the theory begins to become non-perturbative [31] and a small part of the resonance region at m S ∼ m h /2 untested.
Although the relic density and searches at (in)direct detection and collider experiments only probe the mass m S and the portal coupling λ hS , the quartic self-coupling λ S of the scalar singlet does become relevant in the context of DM self-interactions (e.g. [57]) and 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.

Z 3 -symmetric model
The symmetry group that stabilises S is not necessarily Z 2 . We will also consider a complex scalar singlet charged under a Z 3 symmetry, with S transforming as S → e 2π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 [46,58,59], i.e. processes where two DM particles annihilate to an SM particle and another DM particle. 1 The most general scalar potential respecting the Z 3 and SM symmetries is given by 1 It is also possible to have an S 3 term if the Lagrangian is not symmetric under any Zn symmetry. However, such a model also requires quite some tuning to keep the DM sufficiently metastable so that its lifetime is long compared to the age of the Universe. where S † denotes the Hermitian conjugate of S. Unlike the 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 Z 2 -symmetric theory, but has been studied in the context of neutrino masses [60] and in terms of DM phenomenology [6]. The latter included constraints from vacuum stability and perturbativity along with the relic density, direct detection and invisible Higgs decays. Singlet masses below ∼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 Z 2 model is excluded. However, as we will discuss in Section 4.2, vacuum stability sets a limit on µ 3 and thus on the strength of semi-annihilations, so eventually this model also comes within reach of tonne-scale direct detection experiments.

Parameters and nuisances
In Ref.
[44], we studied the direct phenomenological implications of the 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.
[44] 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 gives the parameters of the scalar singlet models and the priors on them that we 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 GeV to 10 TeV, and a second focussed on masses at and below the Higgs resonance m S ∼ 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 M S mass m S (at scale m 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 m 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 crosssections. We thus demand perturbativity as part of the likelihood analysis, by invalidating points in parameter space where any of the dimensionless coupling parameters exceed √ 4π. The choice of this value is similar to that in other studies [52,61].
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 M S mass m h = −2µ 2 H , defined at the scale m S . We then compute the physical pole mass from the input parameters (see discussion in Sec. 4.1). This relationship is affected by radiative corrections from the scalar singlet mass, so the relationship between m h and the pole mass is not constant throughout the parameter space, and we must therefore scan a large range for m h . The resultant value for the pole mass m h is constrained by the likelihood function described in Sec. 4.7.
We scan over a range of ±3σ 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 mean , and the Galactic escape velocity at the solar position v esc . We use the same parameter to control v mean and the rotation speed, v 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 ρ 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 Sec. 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 [31]. In Ref.
[44], 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 mean and the local Galactic escape speed v 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 Z 2 and Z 3 models.
The reduction in the total number of nuisance parameters here compared to Ref.
[44] 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.
[44] with a much less constrained M S 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.
[44]. Our adopted masses for the Z boson, τ lepton and the other quarks, as well as our chosen Fermi and electromagnetic couplings, come from the 2017 compilation of the Particle Data Group [62] (Table 3).

Scanning procedure
Although many directions in parameter space are well constrained, efficient sampling of both the Z 2 and Z 3 scalar singlet models still requires sophisticated sampling algorithms. We scan the parameter space with a differential evolution sampler Diver [23], and an ensemble Markov Chain Monte Carlo (MCMC) T-Walk [23]. 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 S ≈ 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 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 S ∈ [61.8, 63.1] GeV.
We repeat all scans of the Z 3 model with both flat and log priors on µ 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 (λ hS 0.2). We also carry out additional scans with Diver using the older 2017 XENON1T constraint [47] instead of the recent 2018 result [48], 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 [23]. The Diver scans that we present here each used 3400 Intel Xeon Phi 7250 (Knights Landing) cores, for approximately 86 hr in total across all scans. For the 6 T-Walk scans, instead of using a fixed tolerance associated with the sqrtR 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 hr, using the timeout_mins 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 µ 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 × 10 7 and 1.9 × 10 8 samples for the Z 2 and Z 3 models, respectively. With the requirement of absolute vacuum stability, the profile likelihoods are based on 3.3×10 7 and 6.4×10 7 samples for the Z 2 and Z 3 models, respectively.
We produce posteriors and profile likelihoods with pippi [63], basing our posteriors on the maximum posterior density requirement.
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. [67] and refined in Ref. [14]. This implements a matching and running procedure for effective field theories, which is appropriate when m S m t , while not compromising the precision of the Higgs pole mass calculation (due to normal EFT uncertainties from missing O(p 2 /m 2 S ) terms) when m 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 m S . At m S we perform a matching between the SM and the scalar singlet DM model 2 using the FlexibleEFTHiggs matching conditions given in Ref. [14]. At scales larger than m S , we use two-loop RGEs for the scalar singlet model.
For the Z 2 model, the inputs to the FlexibleSUSY spectrum generator are the M S Lagrangian parameters λ hS , λ S , µ H and µ S , defined at the renormalisation scale Q = m S . For the Z 3 -invariant version, this parameter set is extended to include µ 3 (m S ). The parameters λ hS (m S ), λ S (m S ) and µ 3 (m S ) are obtained directly from the model parameters sampled by Scanner-Bit, while µ H (m S ) and µ S (m S ) are obtained by inverting m 2 h = −2µ 2 H and Eq. 2 respectively. 3 FlexibleSUSY fixes the remaining M S parameter λ 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) [76]. These are given in Tables 2 and 3.
Once FlexibleSUSY has determined a consistent set of M S parameters, it computes the singlet and Higgs 2 When mS < mt, we instead set the matching scale to mt, to avoid matching to the UV theory below the scale where we extract SM parameters. 3 In the inversion of Eq. 2 we approximate the M S VEV as Table 3. This means we are effectively making a very mild approximation in the prior for mS, but no such approximation is made in the spectrum calculation, and the impact on the result is negligible.
pole masses, using and Eq. 3. Here Σ S and Σ H incorporate the full oneloop self-energies generated with the help of SARAH 4.12.2 [15][16][17][18]. Σ H additionally includes all two-loop contributions from dominant orders, O(α 2 t +α t α 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 M S mass m h . As the value of the scalar singlet mass can have a significant impact on the relationship between m h (m S ) and m h , we allow m h to vary from 80-180 GeV. This is sufficient to permit a value of 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 ± 0.24 GeV; see Sec. 4.7).
Note that in this setup neither the Higgs pole mass m h , nor the quartic coupling λ h are inputs. Nonetheless, m h is well constrained by the likelihood, so it is important that we can calculate both m h and λ h consistently using the procedure described above.

Vacuum stability and perturbativity
With the running M S 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 Pl = 1.22 × 10 19 GeV and test for vacuum stability.
We classify the stability of the electroweak vacuum in three possible ways: 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 λ(Q 0 ) < 0 for any Q 0 < M Pl , but the electroweak vacuum has an expected lifetime that exceeds the age of the Universe. Unstable. If λ(Q 0 ) < 0 for any Q 0 < M 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. [64] and references therein.
At large field values, the potential can be approximated as V ≈ 1 4 λ 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 [77] The rate of bubble nucleation per unit volume per unit time can be estimated from the bounce action using where Λ B , the scale at which λ h is minimised, has been introduced following Ref. [78]. 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 ≈ e 140 /M Pl , and use it to define the volume of the past lightcone, T 4 U . The predicted number of decays in our past lightcone is therefore We can hence define a likelihood contribution for no decay having occurred in our past light cone as 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 where the rate r is given by, For the SM, this gives a predicted lifetime of ∼1.1 × 10 99 years. Because the dominant contribution from a scalar singlet to the running of λ 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. 4 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 Z 3 model, low-scale vacuum stability gives an additional constraint on the µ 3 parameter. If µ 3 is large, the scalar potential can posess 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 µ 3 parameter. We adopt the condition given in Ref. [6] for an absolutely stable SM vacuum, as an upper limit on µ 3 : This constraint can be relaxed slightly by allowing for the possibility of a 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. [6]). 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 µ 2 S < 0 and λ hS is sufficiently large [26]. Due to the nature of the RGEs for the dimensionless scalar couplings, λ hS and λ S , these couplings only grow with scale.
We let Λ P denote the scale where the dimensionless couplings become larger than our upper bound for perturbativity √ 4π ≈ 3.54. If Λ P < m S , we invalidate the point; otherwise, we record the scale Λ 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 (λ hS and λ S ) results in the theory becoming nonperturbative 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 Λ P . Nevertheless, it is important to consider the order in which we apply these constraints when interpreting the results in Sections 5 and 6.

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 [79] dn S dt Here H is the Hubble rate, σv rel is the thermal average of the relative velocity of DM particles times their selfannihilation cross-section, the number density of DM is given by n S , and its equilibrium number density by n S,eq . When semi-annihilation processes are possible, as in the 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 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 → S * h and S * S * → Sh. We can therefore treat S and S * as the same particle in the Boltzmann calculation, by including a factor of 1/2 [6].
where σv rel is the thermally averaged self-annihilation cross-section without semi-annihilations, and σv rel SS→hS is the equivalent for the semi-annihilation channel. We define a semi-annihilation fraction which we record for each sampled point in the Z 3 parameter space. To deal with semi-annihilations in the Z 3 model, we compute Ω S h 2 using micrOMEGAs 3.6.9.2 [22], with the setting fast = true.
As in Ref.
[44], we employ the measured relic density Ω DM h 2 = 0.1188 ± 0.0010 [80] 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 [12] 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 selfconsistently rescale all direct and indirect signals for the thermal S relic density at each point in the parameter space.

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 Dark-Bit, drawing on the DDCalc [19] implementations of the experimental results of LUX [81], PandaX [49,82] and XENON1T [47,48]. 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 p , follows a Poisson distribution where b denotes the expected number of background events within the analysis region. We interpolate between values in pre-calculated tables contained in DD-Calc 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 [19] for each experiment.
In particular, for the recent XENON1T results [48] 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. [48], 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σ 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 PE, 35 PE] and [35 PE,70 PE]. To this end, we convert a given nuclear recoil energy into an expected S1 signal using Fig. 3 of Ref. [48], 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. [83]. Assuming for simplicity that the remaining (subdominant) 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σ) for a non-zero cross-section at large DM masses.

Indirect detection
Searches for anomalous gamma-ray emission in dwarf spheroidal galaxies constrain the DM annihilation crosssection. The expected flux of gamma rays is for an energy bin of width ∆E i ≡ E max,i −E min,i , where σv 0,j ≡ σv j | v→0 ≡ σv j | s→4m 2 S is the partial annihilation cross-section into final state j in the zero-velocity limit, and dN γ,j /dE is the differential photon multiplicity for annihilations into the jth final state.
We use a combination of analytic expressions from Ref. [31] 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, σv 0 , is given by the processes in the top two rows of Fig. 1 for the Z 3 model, and equivalent processes in the Z 2 model with S = S * . The effective cross-sections for annihilation to SM final states in the Z 3 model are a factor of two smaller than in the Z 2 model, accounting for the fact that only particleantiparticle pairs can annihilate. We compute these by scaling the 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 Z 2 model. With the necessary cross-sections computed we then obtain the predicted spectrum dN γ /dE for each model point by using a Monte-Carlo showering simulation, detailed in Ref. [19]. We then use this to compute a combined likelihood for all the dwarf spheroidals in the Fermi-LAT six-year Pass 8 dataset [84]. The details of this likelihood are given in Ref. [44].

Higgs invisible width
When m 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 Z 3 -charged scalar, the decay width of the Higgs to S bosons is where v 0 is the Higgs VEV. In the Z 2 model the final states are identical, so we must include a symmetry factor of 1/2 to avoid double counting, Eqs. (17) and (18) show that constraints on the Higgs invisible width exclude large λ hS for small singlet masses. With SM-like couplings (which the Higgs possesses in the Z 2 and Z 3 models), the upper limit on the invisible branching fraction of the Higgs is 19% at 95% confidence level [86]. We employ the implementation of the full likelihood associated with this result in DecayBit [64].

Additional likelihoods
We also include simple likelihoods for the nuisance parameters varied in our fits (  PrecisionBit [64]. These quantities are well constrained by existing data. We implement a log-normal likelihood for the local DM density, with a central value ofρ 0 = 0.4 GeV/cm 3 (e.g. [87]) and an uncertainty of σ ρ0 = 0.15 GeV cm −3 , where σ ρ0 = ln(1 + σ ρ0 /ρ 0 ). More details can be found in Ref. [19]. We model the speed distribution of DM in the Milky Way as Maxwell-Boltzmann, truncated at the local Galactic escape velocity v 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 rot [88]. We also constrain the escape velocity using a Gaussian likelihood based on v esc = 550 ± 35 km s −1 , derived from measurements of stellar velocities in the RAVE survey [89].
We apply Gaussian likelihoods to the nuclear parameters as well, based on the estimates σ s = 43 ± 8 MeV [90] and σ l = 50 ± 15 MeV [91]. More detailed discussion of our adopted nuclear and velocity likelihoods can be found in Refs. [19,92].
For the Higgs mass, the top quark mass and the strong coupling, we use Gaussian likelihoods based on m h = 125.09±0.24 GeV [62,93], m t = 173.34±0.76 GeV [62,94] and α s (m Z ) = 0.1181 ± 0.0011 [62].  Contour lines indicate 1σ and 2σ confidence regions, and best fit points are indicated with stars. Shading and white contours show the result of including the 2018 XENON1T analysis [48], whereas grey annotations illustrate the impact of using the 2017 analysis [47] instead. Coloured solid lines indicate published limits from PandaX [49] and XENON1T [48], and the dashed line is the projected sensitivity of LZ [85].

No vacuum constraint
In this section, we present global fits of the 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 √ 4π), up to the greater of m S and m t . That is, we demand that Λ P > max(m S , m t ). The profile likelihoods for this scenario are presented in Fig. 2 in the m S -λ hS (top) and m S -λ S (bottom) planes.
The restriction to Λ P > max(m S , m t ) results in a reduction of the volume of the allowed region compared to our results in Ref. [44]. Any model with values of λ hS or λ S greater than √ 4π at the input scale m S violates the perturbativity condition even before RGE running, so our profile likelihoods extend only to √ 4π in λ hS and λ S . In contrast, Ref.
At very large m S and λ hS , 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 (λ hS ) ≈ 0.2, log 10 (m S /GeV) ≈ 3.6 that was previously allowed in Ref. [44]. This consequence of the perturbativity requirement at high m 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 TeV.
The profile likelihood of the scalar quartic coupling λ S is reasonably uniform over the prior range. This is unsurprising, given that λ S has little phenomenological impact in this model; indeed, this is why it was not included in Ref. [44]. 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 σ SI p for the Z 2 scalar singlet, rescaled by the fraction f ≡ Ω S /Ω 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 [44], significant amounts of parameter space are now excluded from the high-mass modes. The new perturbativity constraint removes parameter space at low σ SI p , whereas the advent of XENON1T cuts into the allowed region from above. LZ [85] 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 [48] by comparing the 1 and 2σ 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.
[44] 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σ preferred region at large m S , where the predictions of the model are consistent with the observed excess.

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 λ hS 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 S ∼ m h /2 is almost entirely ruled out except for the very top of the neck region, where a few points are found with λ hS 0.2 and a stable vacuum. This essentially leaves just the  high mass-modes, centred on approximately 100 GeV and 1 TeV, where λ hS is large enough to stabilise the vacuum.
In Fig. 4 we also show the marginalised posterior for the Z 2 model. As in Ref.
[44], 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 S is 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 S and λ hS , 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 λ hS -dependent radiative corrections cause the Higgs pole mass calculation at m 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 √ 4π. 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.
Let us now take a closer look at just how the Z 2 scalar singlet model can satisfy the vacuum stability constraint. As discussed in Sec. 4.2, we apply the vacuum stability condition by excluding points where we can show perturbatively that the electroweak vacuum is metastable, because λ h becomes negative before the Planck scale M Pl = 1.22 × 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 λ h ≥ 0, ii) some couplings simply run to non-perturbative values before M 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 S -λ S parameter plane, choosing the scale by profiling the likelihood over the other parameters (i.e. plotting Λ P for the best-fit points found in the scan at each combination of m S and λ hS ). We plot the value of Λ P only within the 2σ contours, as determined by the profile likelihood. Note that there can exist parameters with a larger value of Λ P that only have a slightly worse L, and are still within 2σ of the best-fit point. As we run the couplings to a maximum scale of 1 × 10 20 GeV (well above the Planck scale, where quantum gravitational effects become important), points with Λ P equal to this value should be interpreted as valid to at least 1 × 10 20 GeV.
A number of observations can be made from Fig. 5. First of all, we note that for λ S 0.7 the scale of perturbativity is always low, as the singlet quartic coupling quickly runs to non-perturbative values. The dependence on m S is more complicated and can be best understood by comparing to Fig. 2. In the low-mass resonance mode (m S ∼ m h /2), λ hS is typically very small and the scale of perturbativity violation can be very high as long as λ S is sufficiently small. The mode at m S ∼ 100 GeV, on the other hand, requires λ hS > 1, which renders the spectrum invalid at scales well below 10 10 GeV irrespective of λ S . In the high-mass region (m S ∼ 1 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 λ hS ).
In Fig. 6 (left), we show Λ P as function of λ hS and m S in the high-mass region, imposing absolute vacuum stability. There is a rough correlation between λ hS and Λ P , which is only broken for the small λ hS tip of the highmass mode, where Λ P decreases rapidly. This is because such small values of λ hS are insufficient to stabilise the electroweak vacuum, so our requirement that λ h not run negative only finds solutions where λ S contributes to the running of λ h . However, since the impact of λ S on the running of λ h is indirect and only mild, large values of λ 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 λ 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 Λ 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 ∼ 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 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 Λ 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 Λ 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 Z 2 model is to stabilise the electroweak vacuum, it must lie in a narrow region with effective nuclear scattering crosssection 5 × 10 −46 < σ SI p · f < 10 −45 cm 2 and mass 600 GeV < m 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 multiton experiments such as LZ [85] will definitively detect or exclude this scenario.

Best-fit point
The best-fit point for our global fit of the Z 2 model is located at λ S = 6.24 × 10 −3 , λ hS = 2.32 × 10 −4 and m S = 62.48 GeV. This point is located in the low-mass resonance region, the electroweak vacuum is metastable with a lifetime of ∼ 1.1 × 10 99 years, the minimum of λ h occurs at ∼3 × 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. [44], and the portal coupling is approximately a factor of three smaller. Given that the profile likelihood is quite flat with respect to λ hS around the best-fit, the difference in λ hS is not significant.  Indeed, we expect this best-fit point to be similar to what we found in Ref.
[44], as the constraint Λ P > max(m S , m t ) and the variation of λ S do not have a significant impact on the phenomenology at small couplings. Nevertheless, we find ∆ ln L = 0.317 relative to the ideal likelihood (where each individual likelihood takes its maximum value), compared to ∆ ln L = 0.107 in Ref. [44]. 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).    Table 6: Individual contributions to the ∆ log-likelihood for the various best-fit points (see Tabs. 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. [12] for more details on the normalisation). The best-fit points are labelled as follows: A represents a fit with the only constraint that Λ P > max(mS, mt), B is a fit with the additional constraint of absolute vacuum stability, C includes the constraint of Λ P > 10 15 GeV and D also includes the requirement that Ω S h 2 be within 1σ of the observed relic density.
When the constraint of absolute vacuum stability is imposed, the location of the best fit necessarily moves away from the resonance region, where λ hS is too small to stabilise the vacuum. In this case we find a bestfit point at λ S = 2.28 × 10 −3 , λ hS = 2.03 and m S = 3.97 TeV. In this case we find ∆ ln 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 Λ P ∼ 26 TeV. This reduces the theoretical appeal of this point.
By excluding all samples with Λ 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 Λ P = 1.0 × 10 15 GeV, at λ S = 6.59 × 10 −4 , λ hS = 7.36 × 10 −1 , m S = 1.93 TeV. This point has ∆ ln L = 1.121, with the largest contributions coming from direct detection likelihoods. This corresponds to a likelihood ratio Λ = 0.448 relative to the overall best-fit point, which places it inside the generally-preferred 1σ 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.  [48], whereas grey annotations illustrate the impact of using the 2017 analysis [47] instead.
Finally, we consider points with relic densities within 1σ of the Planck measured value, with stable vacua, and that remain perturbative to at least 10 15 GeV. The bestfitting of these points is located at λ S = 2.59 × 10 −3 , λ hS = 6.80 × 10 −1 and m S = 1.94 TeV, and has ∆ ln L = 1.455, still within 1σ 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 ∆ ln L as half the "likelihood χ 2 " of Baker & Cousins [95] 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 ≈ 0.4-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 ≈ 0.15-0.3. This decreases to p ≈ 0.1-0.25 when also requiring that the S relic density is within 1σ 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 Z 3 model

No vacuum constraint
We now turn to the Z 3 scalar singlet model. The main difference compared to the Z 2 model is the presence of an additional parameter µ 3 , which has a significant impact on phenomenology because it can lead to semiannihilations. Fig. 8 presents the profile likelihoods in the m S -λ hS parameter plane (left) and in the m S -µ 3 parameter plane (right), based on scans over the full range of m S . Note that the allowed region for µ 3 is constrained by the vacuum stability condition given in Eq. (11), particularly at small singlet masses.
We find that the allowed resonance region in the Z 3 model is practically identical to the corresponding region in the 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 Z 2 and Z 3 models. The allowed parameter region with m S ∼ 200 GeV is substantially larger in the Z 3 model, and extends to much smaller values of λ hS . This difference in shape can be understood by considering the fraction of semi-annihilation. In Fig. 9 we plot the semiannihilation fraction α, defined in eq. (14), within the 2σ confidence regions. As expected, the extended allowed parameter region in the intermediate mass range corresponds to α ≈ 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 λ hS . In other words, the bound Ω S h 2 ≤ Ω DM h 2 can be evaded at much lower values of λ hS , by invoking a large contribution from semi-annihilation. Similarly, at large values of λ hS the contribution from semiannihilation brings the relic density lower than in the Z 2 model, and direct detection constraints are more easily avoided (given that we consistently rescale signals for the local density of singlet particles).  Regions are shown as a function of mS and the spin-independent direct detection crosssection for scattering with protons, rescaled by the predicted relic abundance σ SI p · f , and compared to the exclusion bounds from various direct detection experiments. Contour lines indicate 1σ and 2σ confidence regions, and best fit points are indicated with stars. Shading and white contours show the result of including the 2018 XENON1T analysis [48], whereas grey annotations illustrate the impact of using the 2017 analysis [47] instead. Other lines indicate limits from PandaX [49] and XENON1T [48], and the projected sensitivity of LZ [85].
For larger masses, semi-annihilations become less efficient, as the semi-annihilation fraction is proportional to µ 3 λ 2 hS /m 6 S at leading order [6]. As a result, the shape of the allowed parameter region is similar to the Z 2 model for m S 1 TeV. The likelihood of this region is however much smaller, and is in fact outside the global 1σ confidence region. The reason is that the Z 3 model requires a complex scalar, whereas we have considered a real scalar for the Z 2 model. The coupling λ hS must therefore be a factor of two larger in the 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. 10 illustrates the impact of current and future direct detection experiments on the parameter space of the Z 3 model. Here we have again rescaled σ SI p by the relic density fraction f = Ω s /Ω DM . Similar to the case of the 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 S m h /2. However, the allowed parameter region at m S ∼ 200 GeV, even though it is substantially larger than its counterpart in the Z 2 model, will eventually be fully explored by direct detection searches.

Absolutely stable vacuum
As with the 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 µ 3 , and in the bottom panel for a scan with logarithmic prior on µ 3 . Although the likelihood is maximised for m S ∼ 100 GeV, the majority of the posterior mass is found at higher masses, m S > 1 TeV. The reason is that at larger masses, the relic abundance becomes independent of µ 3 and therefore benefits strongly from marginalisation (rather than profiling) over µ 3 . Comparing the middle and lower panels, this conclusion is independent of the choice of prior on µ 3 . The choice of prior on µ 3 has a relatively small impact overall, essentially just translating into a stronger preference for large m S when taken flat rather than logarithmic, due to the restriction to lower values of µ 3 at lower m S coming from Eq. 11. Fig. 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 µ 3 on the running of λ 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 Λ P is extremely low (less than 10 10 GeV) throughout 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 Z 2 plot (Fig. 4). Z 3 models remaining in this region after absolute vacuum stability is required are just those with the highest values of λ S , allowing them to combine with non-zero values of µ 3 to send λ 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 semiannihilations in this part of parameter space. A large semi-annihilation fraction requires that the coupling µ 3 is larger than about 300 GeV (see Fig. 9). This in turn forces λ 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 Z 3 model are of limited theoretical appeal from the point of view of stabilising the electroweak vacuum.
In the right panel of Fig. 12, we present the profile likelihood after also imposing the requirement Λ > 10 15 GeV. We find that this not only removes the resonance region, 5 but also the intermediate mass range. In fact, only a tiny region around m S ∼ 1 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 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 [48] already excludes the 2σ 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.

Best-fit point
The best-fit point in the Z 3 model with metastability allowed is at λ S = 4.83 × 10 −1 , λ hS = 3.21 × 10 −4 , µ 3 = 11.8 GeV and m S = 62.48 GeV. This point has a lifetime of ∼1.4 × 10 99 years, and a minimum in its Higgs quartic coupling at ∼3 × 10 13 GeV. For this point we find ∆ ln L = 0.318, essentially the same as for the equivalent best-fit point in the Z 2 model. This result is expected, as the semi-annihilation fraction at this point is α = 0.
With the additional constraint of absolute vacuum stability, the best-fit is located at λ S = 3.52, λ hS = 3.50, µ 3 = 537 GeV and m S = 144 GeV. In this case the semi-annihilation fraction is α = 0.72 and we find ∆ ln L = 0.473. Compared to the equivalent point in the Z 2 model, this represents a small improvement due to the contribution from semi-annihilations. However, Λ P at this point is only 168 GeV, due to the large value of λ S , making it less than appealing in a theoretical sense.
Demanding that Λ P ≥ 10 15 GeV, we find a bestfit point with an absolutely stable vacuum, Λ P = 1.29 × 10 15 GeV and α = 0.004. This point is located at λ S = 3.08 × 10 −2 , λ hS = 6.60 × 10 −1 , µ 3 = 342 GeV and m S = 1.31 TeV. This point has ∆ ln L = 3.975, with the dominant contributions coming from the most recent direct detection experiments. This corresponds to a likelihood ratio Λ = 0.026, which places this point more  than 2σ away from the overall best-fit point. In other words, we find considerable tension in the 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σ of the Planck measured value. The best-fit point under these requirements is located at λ S = 1.19 × 10 −1 , λ hS = 5.92 × 10 −1 , µ 3 = 540 GeV and m S = 1.21 TeV. This point has Ω S h 2 = 0.1137 and Λ P = 1.89 × 10 15 GeV as well as ∆ ln 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 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 ≈ 0.4-0.7, while the best fit with vacuum stability has p ≈ 0.3-0.6. Both of these ranges are very similar to those for equivalent constraints on the Z 2 model, despite semi-annihilations opening up a large region of parameter space. For the model with Λ P > 10 15 GeV, we find p ≈ 0.005-0.02, reducing to p ≈ 0.003-0.01 when also imposing the relic density requirement. This illustrates once again that the 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.

Conclusions
In this work we have investigated two realisations of scalar singlet DM: a real scalar stabilised by a Z 2 symmetry and a complex scalar stabilised by a 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 [48]. 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 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 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 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 Z 2 model. However, the relic density constraint is more severe for a complex scalar than for a real scalar, so that the Z 3 model is in fact more tightly constrained. Indeed, the parameter region with m S ∼ 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 highmass 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 Z 3 symmetry), it remains an interesting possibility for real scalars (with a 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.