Global analyses of Higgs portal singlet dark matter models using GAMBIT

We present global analyses of effective Higgs portal dark matter models in the frequentist and Bayesian statistical frameworks. Complementing earlier studies of the scalar Higgs portal, we use GAMBIT to determine the preferred mass and coupling ranges for models with vector, Majorana and Dirac fermion dark matter. We also assess the relative plausibility of all four models using Bayesian model comparison. Our analysis includes up-to-date likelihood functions for the dark matter relic density, invisible Higgs decays, and direct and indirect searches for weakly-interacting dark matter including the latest XENON1T data. We also account for important uncertainties arising from the local density and velocity distribution of dark matter, nuclear matrix elements relevant to direct detection, and Standard Model masses and couplings. In all Higgs portal models, we find parameter regions that can explain all of dark matter and give a good fit to all data. The case of vector dark matter requires the most tuning and is therefore slightly disfavoured from a Bayesian point of view. In the case of fermionic dark matter, we find a strong preference for including a CP-violating phase that allows suppression of constraints from direct detection experiments, with odds in favour of CP violation of the order of 100:1. Finally, we present DDCalc 2.0.0, a tool for calculating direct detection observables and likelihoods for arbitrary non-relativistic effective operators.

While the nature of the DM particles and their interactions remains an open question, it is clear that the viable candidates must lie in theories beyond the Standard Model (BSM). A particularly interesting class of candidates are weakly interacting massive particles (WIMPs) [5]. They appear naturally in many BSM theories, such as supersymmetry (SUSY) [6]. Due to their weak-scale interaction cross-section, they can accurately reproduce the observed DM abundance in the Universe today. So far there is no evidence that DM interacts with ordinary matter in any way except via gravity. However, the generic possibility exists that Standard Model (SM) particles may connect to new particles via the lowest-dimension gaugeinvariant operator of the SM, H † H . It is therefore natural to assume that the standard Higgs boson (or another scalar that mixes with the Higgs) couples to massive DM particles via such a 'Higgs portal' . The discovery of the Higgs boson in 2012 by ATLAS [28] and CMS [29] therefore opens an exciting potential window for probing DM.
Despite being simple extensions of the SM in terms of particle content and interactions, Higgs portal models have a rich phenomenology, and can serve as effective descriptions of more complicated theories . They can produce distinct signals at present and future colliders, DM direct detection experiments or in cosmic ray experiments. In the recent literature, experimental limits on Higgs portal models were considered from Large Hadron Collider (LHC), Circular Electron Positron Collider and Linear Collider searches, LUX and PandaX, supernovae, charged cosmic and gamma rays, Big Bang Nucleosynthesis, and cosmology [36,41,. The lack of such signals to date places stringent constraints on Higgs portal models.
The first global study of the scalar Higgs portal DM model was performed in Ref. [77]. The most recent global fits [78,79] included relic density constraints from Planck, leading direct detection constraints from LUX, XENON1T, Pan-daX and SuperCDMS, upper limits on the gamma-ray flux from DM annihilation in dwarf spheroidal galaxies with the Fermi-LAT, limits on solar DM annihilation from IceCube, and constraints on decays of SM-like Higgs bosons to scalar singlet particles. The most recent [79] also considered the Z 3 symmetric version of the model, and the impact of requiring vacuum stability and perturbativity up to high energy scales.
In this paper, we perform the first global fits of the effective vector, Majorana fermion and Dirac fermion Higgs portal DM models using the GAMBIT package [80]. By employing the latest data from the DM abundance, indirect and direct DM search limits, and the invisible Higgs width, we systematically explore the model parameter space and present both frequentist and Bayesian results. In our fits, we include the most important SM, nuclear physics, and DM halo model nuisance parameters. For the fermion DM models, we present a Bayesian model comparison between the CP-conserving and CP-violating versions of the theory. We also carry out a model comparison between scalar, vector and fermion DM models.
In Sect. 2, we introduce the effective vector and fermion Higgs portal DM models. We describe the constraints that we use in our global fits in Sect. 3, and the details of our parameter scans in Sect. 4. We present likelihood and Bayesian model comparison results respectively in Sects. 5 and 6, and conclude in Sect. 7. Appendix A documents new features included in the latest version of DDCalc. Appendix B contains all the relevant expressions for the DM annihilation rate into SM particles. All GAMBIT input files, samples and best-fit points for this study are publicly available online via Zenodo [81].

Models
We separately consider vector (V μ ), Majorana fermion (χ ) and Dirac fermion (ψ) DM particles that are singlets under the SM gauge group. By imposing an unbroken global Z 2 symmetry, under which all SM fields transform trivially but (V μ , χ, ψ) → −(V μ , χ, ψ), we ensure that our DM candidates are absolutely stable.
Before electroweak symmetry breaking (EWSB), the Lagrangians for the three different scenarios are [51] where L SM is the SM Lagrangian, W μν ≡ ∂ μ V ν − ∂ ν V μ is the vector field strength tensor, λ hV is the dimensionless vector Higgs portal coupling, λ hχ,hψ /Λ χ,ψ are the dimensionful fermionic Higgs portal couplings, and H is the SM Higgs doublet. The fermionic Lagrangians include both CPodd and CP-even Higgs-portal operators, with θ controlling their relative size. The choice cos θ = 1 corresponds to a pure scalar, CP-conserving interaction between the fermionic DM and the SM Higgs field, whereas cos θ = 0 corresponds to a pure pseudoscalar, maximally CP-violating interaction. We discuss a possible ultraviolet (UV) completion of such a model in Sect. 3.7 (see also Refs. [12,23]). Although all operators in the vector DM model have mass dimension four, the model itself is fundamentally nonrenormalisable, as we do not impose a gauge symmetry to forbid e.g. the mass term for the vector field. Processes with large energies compared to the vector DM mass will violate perturbative unitarity: for large momentum, longitudinal modes of the vector propagator become constant and crosssections become divergent. In this study we remain agnostic as to the origin of the vector mass term and the quartic vector self-interaction, however we do consider perturbative unitarity in Sect. 3.7. After EWSB, the Higgs field acquires a non-zero vacuum expectation value (VEV). In the unitary gauge, we can write where h is the physical SM Higgs field and v 0 = ( √ 2G F ) −1/2 = 246.22 GeV is the Higgs VEV. Thus, the H † H terms in Eqs. (1)(2)(3) generate mass and interaction terms for the DM fields. The tree-level physical mass of the vector DM is For the fermion DM models, the pseudoscalar term (proportional to sin θ ) generates a non-mass-type term that is purely quadratic in the DM fields (e.g., ψγ 5 ψ). Therefore after EWSB, to eliminate this term, we perform a chiral rotation of the fermion DM fields through χ → e iγ 5 α/2 χ, ψ → e iγ 5 α/2 ψ , (6) where α is a real, space-time independent parameter. 1 Using the details outlined in the appendix of Ref. [51], we arrive at the following post-EWSB fermion DM Lagrangians where ξ ≡ θ + α, and m χ,ψ = μ χ,ψ + 1 2 In particular, we note that a theory that is CP-conserving before EWSB (cos θ = 1) is still CP-conserving after EWSB (cos ξ = 1). Because the simplest UV completion leads to cos θ = 1, this means the particular choice of cos ξ = 1 is also natural from the UV perspective. 2 In light of this, we compare the viability of a CP-conserving scenario to the most general case with arbitrary ξ in Sect. 6.

Constraints
The free parameters of the Lagrangians are subject to various observational and theoretical constraints. For the case of IceCube 79-string nulike 1.0.6 [94] vector DM, the relevant parameters after EWSB are the vector DM mass m V and the dimensionless coupling λ hV . 3 The post-EWSB fermion Lagrangians contain three free parameters: the fermion DM mass m χ,ψ , the dimensionful coupling λ hχ,hψ /Λ χ,ψ between DM and the Higgs, and the scalarpseudoscalar mixing parameter ξ .
In Table 1, we summarise the various likelihoods used to constrain the model parameters in our global fit. In the following sections, we will discuss both the physics as well as the implementation of each of these constraints.

Thermal relic density
The time evolution of the DM number density n X is governed by the Boltzmann equation [95] dn X dt where n X,eq is the number density at equilibrium, H is the Hubble rate and σ v rel is the thermally averaged crosssection times relative (Møller) velocity, given by where v cms rel is the relative velocity of the DM particles in the centre-of-mass frame, and K 1,2 are modified Bessel functions. In the case of non-self-conjugate DM, the right hand side of Eq. (11) is divided by two.
In the scenarios discussed above, the annihilation process of DM receives contributions from all kinematically accessible final states involving massive SM fields, including neutrinos. Annihilations into SM gauge bosons and fermions are mediated by a Higgs boson in the s-channel; consequently, near the resonance region, where m X m h /2, it is crucial to perform the actual thermal average as defined in Eq. (12) instead of expanding σ v cms rel into partial waves. 4 Moreover, we take into account the important contributions arising from the production of off-shell pairs of gauge bosons W W * and Z Z * [97]. To this end, for 45 GeV ≤ √ s ≤ 300 GeV, we compute the annihilation cross-section into SM gauge bosons and fermions in the narrow-width approximation via where we employ the tabulated Higgs branching ratios (m * h ) as implemented in DecayBit [84]. For fermionic DM, the dimensionful coupling is implied, λ h X ∈ {λ hV , λ hψ /Λ ψ , λ hχ /Λ χ }. The pre-factor P(X ) is given by In particular, we notice that for CP-conserving interactions of a fermionic DM particle, the annihilation cross-section is p-wave suppressed. As shown in Ref. [97], for √ s 300 GeV the Higgs 1-loop self-interaction begins to overestimate the tabulated Higgs boson width in Ref. [98]. Thus, for √ s > 300 GeV (where the off-shell production of gauge boson pairs is irrelevant anyway), we revert to the tree-level expressions for the annihilation processes given in "Appendix B". Moreover, for m X ≥ m h , DM can annihilate into a pair of Higgs bosons, a process which is not included in Eq. (13). We supplement the cross-sections computed from the tabulated DecayBit values with this process for m X ≥ m h . The corresponding analytical expression for the annihilation cross-sections are given in "Appendix B".
Finally, we obtain the relic density of X by numerically solving Eq. (11) at each parameter point, using the routines implemented in DarkSUSY [99,100] via DarkBit.
In the spirit of the EFT framework employed in this work, we do not demand that the particle X constitutes all of the observed DM, i.e., we allow for the possibility of other DM species to contribute to the observed relic density. Concretely, we implement the relic density constraint using a likelihood that is flat for predicted values below the observed one, and based on a Gaussian likelihood following the Planck measured value DM h 2 = 0.1188 ± 0.0010 [4] for predictions that exceed the observed central value. We include a 5% theoretical error on the computed values of the relic density, which we combine in quadrature with the observed error on the Planck measured value. More details on this prescription can be found in Refs. [80,101].
In regions of the model parameter space where the relic abundance of X is less than the observed value, we rescale all predicted direct and indirect detection signals by f rel ≡ X / DM and f 2 rel , respectively. In doing so, we conservatively assume that the remaining DM population does not contribute to signals in these experiments.

Higgs invisible decays
For m X < m h /2, the SM Higgs boson can decay into a pair of DM particles, with rates given by [51] for the vector, Majorana and Dirac DM scenarios, respectively. These processes contribute to the Higgs invisible width inv , which is constrained to be less than 19% of the total width at 2σ C.L. [102], for SM-like Higgs couplings. We take this constraint into account by using the DecayBit implementation of the Higgs invisible width likelihood, which in turn is based on an interpolation of Fig. 8 in Ref. [102]. Beyond the Higgs invisible width, the LHC provides only a mild constraint on Higgs portal models [63].

Indirect detection using gamma rays
Arguably, the most immediate prediction of the thermal freeze-out scenario is that DM particles can annihilate today, most notably in regions of enhanced DM density. In particular, gamma-ray observations of dwarf spheroidal galaxies (dSphs) of the Milky Way are strong and robust probes of any model of thermal DM with unsuppressed annihilation into SM particles. 5 As described in more detail in Ref. [101], the flux of gamma rays in a given energy bin i from a target object labeled by k can be written in the factorised form Φ i · J k , where Φ i encodes all information about the particle physics properties of the DM annihilation process, while J k depends on the spatial distribution of DM in the region of interest. For s-wave annihilation, one obtains Here κ is a phase space factor (equal to 1 for self-conjugate DM and 1/2 for non-self-conjugate DM), (σ v) 0, j is the annihilation cross-section into the final state j in the zerovelocity limit, and d N γ, j /d E is the corresponding differential gamma-ray spectrum. The J -factor in Eq. (19) is defined via a line of sight (l.o.s.) integral over the square of the DM density ρ X towards the target object k, extended over a solid angle Δ k . In our analysis, we include the Pass-8 combined analysis of 15 dwarf galaxies using 6 years of Fermi-LAT data [85], which currently provides the strongest bounds on the annihilation cross-section of DM into final states containing gamma rays. We use the binned likelihoods implemented in Dark-Bit [101], which make use of the gamLike package. Besides the likelihood associated with the gamma-ray observations, given by we also include a term ln L J that parametrises the uncertainties on the J -factors [85,101]. We obtain the overall likelihood by profiling over the J -factors of all 15 dwarf galaxies, as ln L prof.
Let us remark again that for the case of Dirac or Majorana fermion DM with CP-conserving interactions (i.e., ξ = 0), the annihilation cross-section vanishes in the zero-velocity limit. Scenarios with ξ = 0 therefore pay the price of an additional penalty from gamma-ray observations, compared to the CP-conserving case.

Direct detection
Direct searches for DM aim to measure the recoil of a nucleus after it has scattered off a DM particle [106]. Following the notation of Ref. [101], we write the predicted number of signal events in a given experiment as where M is the detector mass, T exp is the exposure time and φ (E) is the detector efficiency function, i.e., the fraction of recoil events with energy E that are observable after applying all cuts from the corresponding analysis. The differential recoil rate d R/d E for scattering with a target isotope T is given by Here ρ 0 is the local DM density, f (v, t) is the DM velocity distribution in the rest frame of the detector, and dσ/dq 2 (q 2 , v) is the differential scattering cross-section with respect to the momentum transfer q = √ 2m T E.
For the vector DM model, the DM-nucleon scattering process is induced by the standard spin-independent (SI) interaction, with a cross-section given by [51] where μ N = m V m N /(m V +m N ) is the DM-nucleon reduced mass and f N is the effective Higgs-nucleon coupling. The latter is related to the quark content of a proton and neutron, and is subject to (mild) uncertainties. In our analysis we treat the relevant nuclear matrix elements as nuisance parameters; this will be discussed in more detail in Sect. 3.6.
In the case of fermionic DM X ∈ {χ, ψ}, the pseudoscalar current Xiγ 5 X induces a non-standard dependence of the differential scattering cross-section on the momentum transfer q (see e.g., Ref. [107]): where A is the mass number of the target isotope of interest, and F 2 (E) is the standard form factor for spin-independent scattering [108]. As the typical momentum transfer in a scattering process is |q| (1 − 100) MeV m X , we note that direct detection constraints will be significantly suppressed for scenarios that are dominated by the pseudoscalar interaction, i.e., for ξ π/2. For both the vector and fermion models, the spin-dependent (SD) cross-section is absent at leading order. Loop corrections are found not to give a relevant contribution to direct detection in the EFT approach, although they may lead to important effects in specific UVcompletions [109][110][111].
For the evaluation of N p in Eq. (22), we assume a Maxwell-Boltzmann velocity distribution in the Galactic rest frame, with a peak velocity v peak and truncated at the local escape velocity v esc . We refer to Ref. [101] for the conversion to the velocity distribution f (v, t) in the detector rest frame. We discuss the likelihoods associated with the uncertainties in the DM velocity distribution in Sect. 3.6.
We use the DarkBit interface to DDCalc 2.0.0 6 to calculate the number of observed events o in the signal regions for each experiment and to evaluate the standard Poisson likelihood where s and b are the respective numbers of expected signal and background events. We model the detector efficiencies and acceptance rates by interpolating between the precomputed tables in DDCalc. We include likelihoods from the new XENON1T 2018 analysis [89], LUX 2016 [86], PandaX 2016 [87] and 2017 [88] analyses, CDMSlite [90], CRESST-II [91], PICO-60 [92], and DarkSide-50 [93]. Details of these implementations, as well as an overview of the new features contained in DDCalc 2.0.0, can be found in "Appendix A".

Capture and annihilation of DM in the Sun
Similar to the process underlying direct detection, DM particles from the local halo can also elastically scatter off nuclei in the Sun and become gravitationally bound. The resulting population of DM particles near the core of the Sun can then induce annihilations into high-energy SM particles that subsequently interact with the matter in the solar core. Of the resulting particles, only neutrinos are able to escape the dense Solar environment. Eventually, these can be detected in neutrino detectors on the Earth [112][113][114].
The capture rate of DM in the Sun is obtained by integrating the differential scattering cross-section dσ/dq 2 over the range of recoil energies resulting in a gravitational capture, as well as over the Sun's volume and the DM velocity distribution. To this end, we employ the newly-developed public code Capt'n General, 7 which computes capture rates in the Sun for spin-independent and spin-dependent interactions with general momentum-and velocity-dependence, using the B16 Standard Solar Model [115] composition and density distribution. We refer to Refs. [116,117] for details on the capture rate calculation. Notice that similar to direct detection, the capture rate is also subject to uncertainties related to the local density and velocity distribution of DM in the Milky Way. As mentioned earlier, these uncertainties are taken into account by separate nuisance likelihoods to be discussed in Sect. 3.6.
Neglecting evaporation (which is well-justified for the DM masses of interest in this study [118][119][120]), the total population of DM in the Sun N X (t) follows from where C(t) is the capture rate of DM in the Sun, and A(t) ∝ σ v rel N X (t) 2 is the annihilation rate of DM inside the Sun; this is calculated by DarkBit. We approximate the thermally averaged DM annihilation cross-section, which enters in the expression for the annihilation rate, by evaluating σ v at v = √ 2T /m X , where T = 1.35 keV is the core temperature of the Sun. 7 https://github.com/aaronvincent/captngen. At sufficiently large t, the solution for N X (t) reaches a steady state and depends only on the capture rate. However, the corresponding time scale τ for reaching equilibrium depends also on σ v, and thus changes from point to point in the parameter space. Hence, we use the full solution of Eq. (27) to determine N X at present times, which in turn determines the normalization of the neutrino flux potentially detectable at Earth. We obtain the flavour and energy distribution of the latter using results from WimpSim [121] included in DarkSUSY [99,100].
Finally, we employ the likelihoods derived from the 79string IceCube search for high-energy neutrinos from DM annihilation in the Sun [94] using nulike [122] via DarkBit; this contains a full unbinned likelihood based on the eventlevel energy and angular information of the candidate events.

Nuisance likelihoods
The constraints discussed in the previous sections often depend on nuisance parameters, i.e. parameters not of direct interest but required as input for other calculations. Examples are nuclear matrix elements related to the DM direct detection process, the distribution of DM in the Milky Way, or SM parameters known only to finite accuracy. It is one of the great virtues of a global fit that such uncertainties can be taken into account in a fully consistent way, namely by introducing new free parameters into the fit and constraining them by new likelihood terms that characterise their uncertainty. We list the nuisance parameters included in our analysis in Table 2, and discuss each of them in more detail in the rest of this section.
Following the default treatment in DarkBit, we include a nuisance likelihood for the local DM density ρ 0 given by a log-normal distribution with central value ρ 0 = 0.40 GeV cm −3 and an error σ ρ 0 = 0.15 GeV cm −3 . To reflect the log-normal distribution, we scan over an asymmetric range for ρ 0 . For more details, see Ref. [101]. For the parameters determining the Maxwell-Boltzmann distribution of the DM velocity in the Milky Way, namely v peak and v esc , we employ simple Gaussian likelihoods. Since v peak is equal to the circular rotation speed v rot at the position of the Sun for an isothermal DM halo, we use the determination of v rot from Ref. [123] to obtain v peak = 240 ± 8 km s −1 . 8 The escape velocity takes a central value of v esc = 533 ± 31.9 km s −1 , where we convert the 90% C.L. interval obtained by the RAVE collaboration [126], assuming that the error is Gaussian.
As noted already in Sect. 3.4, the scattering cross-section of DM with nuclei (which enters both the direct detection and solar capture calculations) depends on the effective DMnucleon coupling f N , which is given by [101] Here f

(N )
T q are the nuclear matrix elements associated with the quark q content of a nucleon N . As described in more detail in Ref. [127], these are obtained from the following observable combinations where m l ≡ (m u + m d )/2. We take into account the uncertainty on these matrix elements via Gaussian likelihoods given by σ s = 43 ± 8 MeV [128] and σ l = 50 ± 15 MeV [129]. The latter deviates from the default choice implemented in DarkBit as it reflects recent lattice results, which point towards smaller values of σ l (see Ref. [129] for more details). Furthermore, we have confirmed that the uncertainties on the light quark masses have a negligible impact on f N . Thus, for simplicity, we do not include them in our fit. We also use a Gaussian likelihood for the Higgs mass, based on the PDG value of m h = 125.09 ± 0.24 GeV [130]. In line with our previous study of scalar singlet DM [78], we allow the Higgs mass to vary by more than 4σ as the phenomenology of our models depends strongly on m h , most notably near the Higgs resonance region. Finally, we take into account the uncertainty on the strong coupling constant α s , which enters the expression for the DM annihilation crosssection into SM quarks (see "Appendix B"), taking a central value α M S s (m Z ) = 0.1181 ± 0.0011 [130].

Perturbative unitarity and EFT validity
The parameter spaces in which we are interested are limited by the requirement of perturbative unitarity. First of all, this requirement imposes a bound on any dimensionless coupling in the theory. Furthermore, as neither the vector or fermion Higgs portal models are renormalisable, we must ensure that the effective description is valid for the parameter regions to be studied. The dimensionless coupling λ hV in the vector DM model is constrained by the requirement that annihilation processes such as V V → hh do not violate perturbative unitarity. Determining the precise bound to be imposed on λ hV is somewhat involved, so we adopt the rather generous requirement λ hV < 10 with the implicit understanding that perturbativity may become an issue already for somewhat smaller couplings.
For small DM masses, an additional complication arises from the fact that theories with massive vector bosons are not generally renormalisable. In that case cross-sections do not generally remain finite in the m V → 0 limit and a significant portion of parameter space violates perturbative unitarity [131]. However, by restricting ourselves to the case of μ 2 V , λ hV ≥ 0 we can safely tackle both issues due to the fact that m V → 0 implies λ hV → 0. Using Eq. (5), this condition translates to A more careful analysis might lead to a slightly larger valid parameter space, but as we will see in Sect. 5.1.1, those regions would be excluded by the Higgs invisible width anyway.
The EFT validity of the fermion DM models depends on the specific UV completion. To estimate the range of validity, we consider a UV completion in which a heavy scalar mediator field Φ couples to the fermion DM X and the Higgs doublet as [12] where X ∈ {χ, ψ} and μ has mass dimension 1. 9 For this specific UV completion, we assume that the mixing between Φ and the Higgs field is negligible and can be ignored. The heavy scalar field can be integrated out to give a dimensionful coupling in the EFT approximation as Thus, by comparing Eq. (32) with the fermion DM Lagrangians in Eqs. (2) and (3), we can identify μg X g H /m 2 Φ with λ h X /Λ X . As μ should be set by the new physics scale, we take it to be roughly m Φ , implying g X g H /m Φ ∼ λ h X /Λ X . In addition, we require the couplings to be perturbative, i.e., g X g H ≤ 4π .
We need to consider the viable scales for which this approximation is valid. We require that the mediator mass m Φ is far greater than the momentum exchange q of the interaction, i.e., m Φ q such that Φ can be integrated out. For DM annihilations, the momentum exchange is q ≈ 2m X . Thus, the EFT approximation breaks down when m Φ < 2m X and our EFT assumption is violated when As the typical momentum transfer in a direct detection experiment is roughly on the order of a few MeVs, the EFT validity limit requires m Φ O (MeV), which is always satisfied by the previous demand m Φ > 2m X for the mass ranges of interest. In this case, we assume that the couplings saturate the bound from perturbativity, i.e., g X g H = 4π ; the constraint would be stronger if the couplings were weaker.
For parameter points close to the EFT validity bound, the scale of new physics is expected to be close to or even below 2m χ . In this case, the annihilation cross-section σ v rel , used in predictions of both the relic density and indirect detection signals, may receive substantial corrections from interactions with Φ, which are not captured in the EFT approach. The likelihoods computed for these points should hence be interpreted with care.
Note that this prescription is only the simplest and most conservative approach; additional constraints can be obtained by unitarising the theory (e.g. [132]).

Scan details
We investigate the Higgs portal models using both Bayesian and frequentist statistics. The parameter ranges and priors that we employ in our scans of the vector and fermion DM models are summarised in Tables 3 and 4, respectively. Whilst the likelihoods described in the previous sections are a common ingredient in both our frequentist and Bayesian analyses, the priors only directly impact our Bayesian analyses. We discuss our choice of priors in Sect. 5.2. For a review of our statistical approaches to parameter inference, see e.g., Ref. [80].
There are two main objectives for the Bayesian scans: firstly, producing marginal posteriors for the parameters of interest, where we integrate over all unplotted parameters, and secondly, computing the marginal likelihood  Table 4 Parameter ranges and priors for the fermion DM models. Our choice for the range of ξ between 0 and π reflects the fact that only odd powers of cos ξ appear in the observables that we consider, but never odd powers of sin ξ , which cancel exactly due to the complex conjugation. Thus, the underlying physics is symmetric under ξ → −ξ For the frequentist analysis, we are interested in mapping out the highest likelihood regions of our parameter space. For this analysis we largely use Diver, a differential evolution sampler, efficient for finding and exploring the maxima of a multi-dimensional function. Details of T-Walk and Diver can be found in Ref. [136].
Due to the resonant enhancement of the DM annihilation rate by s-channel Higgs exchange at m X ≈ m h /2, there is a large range of allowed DM-Higgs couplings that do not overproduce the observed DM abundance. When scanning over the full mass range, it is difficult to sample this resonance region well, especially with a large number of nuisance parameters. For this reason, we perform separate, specific scans in the low-mass region around the resonance, using both T-Walk and Diver. When plotting the profile likelihoods, we combine the samples from the low-and high-mass scans.
In addition, as part of the Bayesian analysis, we perform targeted T-Walk and MultiNest scans of the fermion DM parameter space where the interaction is wholly scalar (ξ = 0), using the same priors for the fermion DM mass and its dimensionful coupling as in Table 4. This allows us to perform model comparison between the cases where ξ is fixed at zero, and where it is left as a free parameter. The convergence criteria that we employ for the different samplers are outlined in Table 5. We carried out all Diver scans on 340 Intel Xeon Phi 7250 (Knights Landing) cores. As in our recent study of scalar singlet DM [79], we ran T-Walk scans on 1360 cores for 23 h, providing us with reliable sampling. The MultiNest scans are based on runs using 240 Intel Broadwell cores, with a relatively high tolerance value, which is nevertheless sufficient to compute the marginal likelihood to the accuracy required for model comparison. We use the importance sampling log-evidence from MultiNest to compute Bayes factors.
For profile likelihood plots, we combine the samples from all Diver and T-Walk scans, for each model. The plots are based on 1.46 × 10 7 , 1.70 × 10 7 and 1.73 × 10 7 samples for the vector, Majorana and Dirac models, respectively. We do all marginalisation, profiling and plotting with pippi [137].

Profile likelihoods
In this section, we present profile likelihoods from the combination of all Diver and T-Walk scans for the vector, Majorana and Dirac models. Profile likelihoods in the vector model parameters are shown in Fig. 1, with key observables rescaled to the predicted DM relic abundance in Fig. 2. Majorana model parameter profile likelihoods are shown in Figs. 3 and 4, with observables in Fig. 5. For the Dirac model, we simply show the mass-coupling plane in Fig. 6, as the relevant physics and results are virtually identical to the Majorana case. Figure 1 shows that the resonance region is tightly constrained by the Higgs invisible width from the upper-left when m V < m h /2, by the relic density constraint from below, and by direct and indirect detection from the right. Nevertheless, the resonant enhancement of the DM annihilation at around m h /2, combined with the fact that we allow for scenarios where V μ is only a fraction of the observed DM, permits a wide range of portal couplings. Interestingly, the perturbative unitarity constraint (shown as dark grey) in Eq. (30) significantly shortens the degenerate 'neck' region that appears exactly at m h /2. Most notably, this is in contrast with the scalar Higgs portal model [78,79] where no such constraint exists.

Vector model
The high-mass region contains a set of solutions at m V 10 TeV and λ hV 1, which are constrained by the relic density from below and direct detection from the left. This  second island is prominent in both our previous studies of the scalar Higgs portal model [78,79] as well as other studies of the vector Higgs portal [51]. The precise extent of this region depends on the upper bound imposed on λ hV to reflect the breakdown of perturbativity. While the constraint that we apply ensures that perturbative unitarity is not violated [131], higher-order corrections may nevertheless become important in this region. The perturbative unitarity constraint from Eq. (30) excludes solutions that would otherwise exist in a separate triangular region at m χ m h , λ hV 1. In Table 6, we show a breakdown of the contributions to the likelihood at the best-fit point, which lies on the lower end of the resonance region at λ hV = 4.9 × 10 −4 and m V = 62.46 GeV. If we demand that vector singlet DM constitutes all of the observed DM, by requiring V h 2 to be within 1σ of the observed Planck relic abundance, the best-fit point remains almost unchanged, at λ hV = 4.5 × 10 −4 and m V = 62.46 GeV. We give details of these best-fit points, along with the equivalent for fermion models, in Table 7.
In Fig. 2, we show the relic density of the vector model (top), as well as the cross-sections relevant for direct (centre) and indirect detection (bottom), all plotted as a function of mass. Only models along the lower-λ hV edge of the two likelihood modes have relic densities equal to the observed value. Larger values of λ hV result in progressively larger annihilation cross-sections and therefore more suppression of the relic density, cancelling the corresponding increase in σ SI p and resulting in an essentially constant rescaled crosssection f ·σ SI p ∼ 10 −45 cm −2 in the remaining allowed highmass region. Future direct detection experiments such as LZ [138] will be able to probe the high-mass region in its entirety. However, the best-fit point -near the bottom of the resonance region -will remain out of reach. Future indirect searches, such as the Cherenkov Telescope Array (CTA) 10 [139] will also be able to probe large amounts of the high-mass region; however it does not have the exclusion power that direct detection does for Higgs portal models. Again, the best-fit point remains out of reach.

Majorana fermion model
We show profile likelihoods in the (m χ , λ hχ /Λ χ ) plane in Fig. 3, with the low-mass region in the left panel and the full mass region in the right panel. Here, there are no longer two distinct solutions: the resonance and high mass regions are connected. From the left panel in Fig. 4, where we plot the profile likelihood in the (m χ , ξ) plane, we can see that these regions are connected by the case where the portal interac- 10 The CTA projections plotted in Fig. 2 assume an Einasto density profile, and are based on 500 h of observations of the Galactic centre [139], with no systematic uncertainties. They should therefore be considered optimistic [140,141].   The high mass region prefers ξ ∼ π/2, with a wider deviation from π/2 permitted as m χ is increased, due to direct detection constraints, which become less constraining at higher WIMP masses. There is an enhancement in the permitted range of mixing angles at m χ m h , due to the contact term (∝ χχhh), where DM annihilation to on-shell Higgses reduces the relic density, providing another mechanism for suppressing direct detection signals, thus lifting the need to tune ξ .
The results are roughly symmetric about ξ = π/2, however due to odd powers of cos ξ in the annihilation cross-section (see "Appendix B"), there is a slight asymmetry for masses above m h . This is most clearly seen in the triangular 'wings' at m χ m h in Fig. 4 where there are more solutions for ξ > π/2 than for ξ < π/2.
In the resonance region, we see the same triangular region as in the vector DM case: bounded from below by the relic density, and from the upper-left by the Higgs invisible width. However, in contrast to the vector DM case where direct detection limits squeeze the allowed region from the upper right, the addition of the mixing angle ξ as a free parameter allows for the fermionic DM models to escape these constraints. As the pseudoscalar coupling is increased and the scalar coupling is correspondingly decreased, the SI crosssection becomes steadily more q 2 -suppressed, as seen in  Eq. (25). Noting that, the neck region at m χ = m h /2 is less well-defined than in the vector and scalar DM cases above the triangle region. Notably however, as the SI cross-section becomes steadily more q 2 -suppressed, the annihilation crosssection becomes less p-wave suppressed (Eq. 14), and indirect detection comes to dominate the constraint at the edge of the allowed parameter space just above the resonance.
In the low-mass resonance region, virtually all values of the mixing angle are permitted, seen clearly in the left panel of Fig. 4, as even purely scalar couplings are not sufficient for direct detection to probe the remaining parameter space. The right panel also shows this in the lower 'bulb': couplings between 10 −3 and 10 −5 GeV −1 are only permitted in the resonance region, without any constraint on the mixing angle.
In the high-mass region, we see that unlike the vector DM case, a wide range of WIMP masses between 100 GeV and 10 TeV are acceptable, with degenerate maximum likelihood. This is again due to the q 2 -suppression of the direct detection constraints when considering all possible values of ξ . The large triangular high-mass region is constrained by the EFT validity constraint from above (highlighted in dark grey) and the relic density constraint from below.
In Fig. 5, we show the relic density (top) and scaled cross-sections for direct (centre) and indirect detection (bottom). For plotting purposes, we compute σ SI at a reference momentum exchange of q = 50 MeV, typical of direct detection experiments. Substantial fractions of allowed parameter space lie close to current limits, but unsurprisingly, large portions of the parameter space will not be probed by future direct detection experiments, due to the momentum suppression. This is also true for indirect detection, where crosssections are velocity suppressed. However, given that the two suppressions have opposite dependences on the mixing parameter, the two probes will be able to compensate for each others' weaknesses to a certain extent. Table 6 shows a breakdown of the contributions to the likelihood at the best-fit point, which lies in the high mass region at m χ = 138.4 GeV, λ hχ /Λ χ = 4.5 × 10 −2 GeV −1 and ξ = 1.96 rad (Table 7). When we demand that χ saturates the observed DM relic abundance, the best-fit point shifts to the lower end of the resonance region at m χ = 61.03 GeV, λ hχ /Λ χ = 6.3 × 10 −6 GeV −1 and ξ = 1.41 rad.

Dirac fermion model
The results from our low-and high-mass scans of the Dirac fermion model are very similar to those for the Majorana model. We therefore only show results in the (m ψ , λ hψ /Λ ψ ) plane in Fig. 6.
In Table 6  Log-likelihood contribution up all of the DM, the best-fit point shifts slightly to the bottom of the high mass triangle at λ hψ /Λ ψ = 3.6 × 10 −4 GeV −1 , m ψ = 9.9 TeV and ξ = 2.07 rad. We compare the locations of these best-fit points to those from the vector and Majorana models in Table 7.

Goodness of fit
In Table 6, we show the contribution to the log-likelihood for the best-fit points of the vector, Majorana and Dirac DM models. By equating Δ ln L to half the "likelihood χ 2 " of Baker and Cousins [142], we can compute an approximate p-value for each best-fit point against a null hypothesis. We take this null to be the 'ideal' case, which we define as the background-only contribution in the case of exclusions, and the observed value in the case of detections.
For the vector DM model, using either one or two effective degrees of freedom, we find a p-value between roughly 0.4 and 0.7. Requiring the relic density of V μ to be within 1σ of the Planck value, the p-value becomes p ≈ 0.35-0.65. For both the Majorana and Dirac fermion models, we also find p ≈ 0.4-0.7, falling to 0.35-0.65 with the relic density requirement. All of these are completely acceptable p-values.

Marginal posteriors
The marginal posterior automatically penalises fine-tuning, as upon integration of the posterior, regions with a limited 'volume of support' over the parameters that were integrated over are suppressed. 11 As usual, the marginal posteriors depend upon the choice of priors for the free model parameters, which are summarised in Tables 3 and 4. We choose flat priors where parameters are strongly restricted to a particular scale, such as the mixing parameter and the DM mass in scans restricted to the low-mass region. For other parameters, in order to avoid favouring a particular scale we employ logarithmic priors. Note that in this treatment for the fermionic DM models we have not chosen priors that favour the CP-conserving case. We instead present posteriors for this well motivated case separately, and later in Sect. 6 we perform a Bayesian model comparison between a CP-conserving fermionic DM model and the full model considered here.

Vector model
To obtain the marginal posterior distributions, we perform separate T-Walk scans for the low and high mass regimes, 11 By 'volume of support', we mean the regions of the parameter space that have a non-negligible likelihood times prior density. shown in Fig. 7. Within each region we plot the relative posterior probability across the parameter ranges of interest.
In the left panel of Fig. 7, the scan of the resonance region shows that the neck region is disfavoured after marginalising over the nuisance parameters, particularly m h , which sets the width of the neck. This dilutes the allowed region due to volume effects.
In the full-mass-range scan, the fine-tuned nature of the resonance region is clearly evident. Although the best-fit point in the profile likelihood lies in the resonance region, the posterior mass is so small in the entire resonance region that it drops out of the global 2σ credible interval.

Majorana fermion model
As already seen in the profile likelihoods, for the case of Majorana fermion DM, the presence of the mixing parameter ξ leads to a substantial increase in the preferred parameter region (see Fig. 8). In the resonance region (left panel), there is now a thin neck-like region at m χ ≈ m h /2. This neck region is the same one seen in both the scalar and vector profile likelihoods, but falls within the 2σ credible region of the Majorana posterior, as the admittance of ξ reduces direct detection constraints (Eq. 25), softening the penalisation from integrating over nuisance parameters. When we compute the posterior over the full mass range, we once again find the resonance region to be somewhat disfavoured, but now there are large parameter regions with high posterior probabilities for m χ > m h .
Nevertheless, direct detection does have a significant impact on the high-mass region, in spite of the mixing parameter ξ . While the 2σ contour is roughly triangular, the points with highest posterior probability (i.e. within the 1σ contours) are split into two smaller triangles. The approximately rectangular region that separates these two triangular regions is disfavoured by the combination of volume effects and direct detection, which requires ξ to be tuned relatively close to π/2.  To better understand the role of tuning in ξ in the process of marginalisation, we show the marginalised posterior in the (m χ , ξ) and (ξ, λ hχ /Λ χ ) planes in Figs. 9 and 10, respectively. Figure 9 provides a clear understanding of the differences between the marginalised posteriors in Fig. 8 and the profile likelihood in Fig. 3. In the resonance region (left panel), the neck region is less prominent in the marginalised posterior because direct detection limits become very constraining as soon as m χ > m h /2 and the mixing parameter is forced to be very close to π/2. In the full-range scan (right panel) we see the annihilation channel χχ → hh open up, thus allowing a greater range of values for ξ , leading to an enhancement in the marginalised posterior probability. This clearly corresponds to the 1σ triangular region in the masscoupling plane at m χ ≈ m h , in the right hand panel of Fig. 8.
In the left panel of Fig. 10, which focuses on the resonance region, we see two separate solutions for the mixing angle and coupling: the larger island at lower coupling corresponds to the triangular region at m χ < m h /2, permitting all values of ξ , and the thinner solution at larger couplings reflects the solution at m χ > m h /2, where the scalar coupling between the Higgs and the Majorana DM needs to be suffi-  ciently small (i.e. ξ ∼ π/2) to evade direct detection limits. The two regions appear disconnected because the intermediate parameter points require so much tuning that they fall outside of the 2σ credible regions upon marginalisation. Considering the full mass range (see the right panel in Fig. 10), we find that the lower 'bulb' seen in the profile likelihood in Fig. 4 is hardly visible in the marginalised posterior when integrating over the nuisance parameters, due to a lower posterior volume in the resonance region. We can condense the information from Figs. 9 and 10 further by marginalising over all parameters except for ξ , thus obtaining a 1D posterior probability. The result is shown in Fig. 11, where the preference for ξ ≈ π/2 becomes clear. In other words, for the case of Majorana fermion DM, there is a strong preference for permitting an increased admixture of pseudoscalar-type couplings to suppress the constraints from direct detection and the relic density, due to a momentum and velocity suppressed cross-section respectively.
For comparison, we consider the CP-conserving case with fixed ξ = 0 in Fig. 12. As expected from the discussion above, we find that the permitted parameter space shrinks vastly with respect to the case where the mixing parameter is allowed to vary (see Fig. 8). In the resonance region (left panel), we see that direct detection, the invisible Higgs width and relic density impose strong constraints from the left, upper-left and below, respectively. No neck region exists because the direct detection constraints are too strong, overlapping with constraints on the invisible width of the Higgs boson. In the full-range scan (right panel), we find that the only surviving parameter space is split into the resonance These islands are constrained by direct detection and the EFT validity requirement. Both will be ruled out by the next generation of direct detection experiments, if no DM signal is observed. Our analysis of the Dirac fermion model parameter space is identical to the Majorana fermion one, whether ξ is fixed or left as a free parameter, so to avoid repetition we omit those results.
It should be clear from the comparison between Figs. 8 and 12 that the CP-conserving case (ξ = 0) is strongly disfavoured relative to the case where ξ is allowed to vary.
We will make this qualitative observation more precise in the following section.

Background
To be able to comment on the relative plausibility of the different Higgs portal models, we must also perform a quantitative model comparison. To do this, we compute Bayes factors for pairs of models, say M 1 and M 2 as [143][144][145], where Z(M) is the evidence of a model M. This is the integral of the likelihood of the observed data L(D|θ) over the possible parameter values θ in that model, weighted by the prior on the parameters P(θ ), We perform this integration using MultiNest [133,134], which is designed to calculate the Bayesian evidence. The final odds ratio (of the probability that M 1 is correct to the probability that M 2 is correct) is the product of the Bayes factor and the ratio of any prior beliefs P(M 1 )/P(M 2 ) that we might have in these models,  In our analysis, we take the prior probability of every model to be equal such that the factor, for all pairs of models. Thus, the odds ratio is simply given by the Bayes factor. Note that even when computing the evidence for or against CP violation in the fermionic model below, we do not apply any further prior in favour of CP conservation. The volume integrals involved in the Bayes factor automatically implement the concept of naturalness via Occam's razor, penalising models with more free parameters if they do not fit the observed data any better than models with less parameters. From Eq. (35), we can see that the evidence of a model depends on the prior choices for its parameters. This prior on the model parameters (along with the priors on the models themselves) makes the results of Bayesian model comparison inherently prior-dependent. However, the influence of common parameters treated with identical priors in both models approximately cancels when taking the ratio of evidences, as in Eq. (34). The overall prior dependence of the Bayes factor can thus be minimised by minimising the number of non-shared parameters between the models being compared. The best case is where one model is nested inside the other, and corresponds simply to a specific choice for one of the degrees of freedom in the larger model. In this case, the leading prior dependence is the one coming from the chosen prior on the non-shared degree of freedom. Thus, we first investigate the question of CP violation in the Higgs portal, which we can address in this manner, before going on to the more prior-dependent comparison of the broader models.

CP violation in the Higgs portal
We perform Bayesian model comparison for the fermionic Higgs portal DM, and nested variants of it, by comparing the CP-conserving case (ξ = 0) to the model where the CP phase of the portal coupling is allowed to vary freely. Due to the similarity of the likelihood for the Dirac and Majorana fermion models, we do this for the Majorana fermion model only. We carry out this exercise for two different parametrisations of the model, corresponding to two different priors on the larger parameter space in which the CP-conserving scenario is nested: 1. Assuming the parametrisation that we have discussed thus far for the Majorana model, taking a uniform prior for ξ and a logarithmic prior for λ hχ /Λ χ . This corresponds to the assumption that some single mechanism uniquely determines the magnitude and phase of both couplings.
2. Assuming that the scalar and pseudoscalar couplings originate from distinct physical mechanisms at unrelated scales, such that they can be described by independent logarithmic priors. The post-EWSB Lagrangian in this parametrisation contains the terms In this case, the parameters ξ and λ hχ /Λ χ from the first parametrisation are replaced by g s /Λ s and g p /Λ p . In this parametrisation, the Bayes factor may be sensitive to the range of the prior for the couplings, as the normalisation factor does not cancel when computing the Bayes factor for the CP-conserving scenario. We choose −6 ≤ log 10 (g/Λ) ≤ 0 for the couplings when computing the Bayes factors in this parametrisation, in line with the prior that we adopt for λ hχ /Λ χ in parametrisation 1.
The CP-conserving model is nested within both of these models, as ξ = 0 in the first, and as g p /Λ p = 0 in the second (although the exact limit of ξ = 0 is not contained within our chosen prior for the second parameterisation, seeing as we choose a logarithmic prior on g p ). As stated in Eq. 37, the ratio of the prior probabilities for these models is taken to be 1 here, and is not related to priors of parameters discussed above. We are comparing two separate models: one with a pure CP-even coupling between the DM fermion and the Higgs and another model where there is also a pseudoscalar coupling, which a priori is very unlikely to be zero.
In Table 8, we give the odds ratios against the CPconserving case in each of these parametrisations. The value given in the final column of this table is the ratio of the evidence for the CP-violating model to the CP-conserving case. Depending on the choice of parametrisation, we see that there is between 140:1 and 70:1 odds against the CP-conserving version of the Majorana Higgs portal model. The similarity in order of magnitude 12 between these two results is expected, as it reflects the relatively mild prior-dependence of the Bayes factor when performing an analysis of nested models that differ by only a single parameter. Given the similarity of the likelihood functions in the Majorana and Dirac fermion models, the odds against the pure CP-conserving version of the Dirac fermion Higgs portal model can also be expected to be very similar.
The odds ratio tells us the relative plausibility of one model relative to the other. According to the standard scale frequently used for interpreting Bayesian odds ratios (the Jeffreys scale; [143,144]), this constitutes strong evidence against pure CP-even coupling in fermionic Higgs portal models. The preference for a CP-violating coupling can also be seen in Fig. 11, where there is a clear preference for ξ = π/2, whereas the CP-even coupling falls outside of the 2σ credible region.

Scalar, vector, Majorana or Dirac?
We also carry out model comparison between the different Higgs portal models: scalar, vector, Majorana and Dirac. As these models are not nested, they each have unique parameters. This means that there is no a priori relationship between their respective parameters that would allow the definition of equivalent priors on, e.g., masses or couplings in two different models. The prior dependence of the Bayes factor is therefore unsuppressed by any approximate cancellations when taking the ratio of evidences in Eq. (34). We caution that the resulting conclusions are consequently less robust than for the nested Majorana models. For this exercise, we update the fit to the scalar model from Ref. [78] to incorporate the likelihood function and nuisances that we use in the current paper. We find that the scalar Higgs portal model has the largest evidence value in our scans, but is comparable to the fermion DM models. In Table 9, we give the odds ratios against each of the Higgs portal models, relative to the scalar model. The data have no preference between scalar and either form of fermionic Higgs portal model, with odds ratios of 1:1. The vector DM model is disfavoured with a ratio of 6:1 compared to the scalar and fermion models; this constitutes 'positive' evidence against the vector DM model according to the Jeffreys scale, though the preference is only rather mild. Overall, there is no strong preference for Higgs portal DM to transform as a scalar, vector or fermion under the Lorentz group.
As we find no strong preference between the different Higgs portal DM models using logarithmic priors, we omit a dedicated prior sensitivity analysis. If different assumptions on priors were to yield a stronger preference for any of the models under consideration, the only conclusion would be that such a preference is not robust to changes in the prior. The situation is hence different from the one in Sect. 6.2, where we did find a strong preference against the CP-conserving model, which we showed to be largely independent of the assumed prior.

Conclusions
In this study we have considered and compared simple extensions of the SM with fermionic and vector DM particles stabilised by a Z 2 symmetry. These models are nonrenormalisable, and the effective Higgs-portal coupling is the lowest-dimension operator connecting DM to SM particles. Scenarios of this type are constrained by the DM relic density predicted by the thermal freeze-out mechanism, invisible Higgs decays, and direct and indirect DM searches. Perturbative unitarity and validity of the corresponding EFT must also be considered. We find that the vector, Majorana and Dirac models are all phenomenologically acceptable, regardless of whether or not the DM candidate saturates the observed DM abundance. In particular, the resonance region (where the DM particle mass is approximately half the SM Higgs mass) is consistent with all experimental constraints and challenging to probe even with projected future experiments. On the other hand, larger DM masses are typically tightly constrained by a combination of direct detection constraints, the relic density requirement and theoretical considerations such as perturbative unitarity. Our results show that with the next generation of direct detection experiments (e.g., LZ [138]), it will be possible to fully probe the high-mass region for both the vector and CPconserving fermion DM model. Future indirect experiments such as CTA [139] will be sensitive to parts of viable parameter space at large DM masses, but will have difficulty in probing the resonance region.
An interesting alternative is fermionic DM with a CPviolating Higgs portal coupling, for which the scattering rates in direct detection experiments are momentum-suppressed. By performing a Bayesian model comparison, we find that data strongly prefers the model with CP violation over the CP-conserving one, with odds of order 100:1 (over several priors). This illustrates how increasingly tight experimental constraints on weakly-interacting DM models are forcing (2019) 79:38 Page 21 of 28 38 us to abandon the simplest and most theoretically appealing models, in favour of more complex models. We have also used Bayesian model comparison to determine the viability of the scalar Higgs portal model relative to the fermionic and vector DM models. We find a mild preference for scalar DM over vector DM, but no particular preference between the scalar and the fermionic model. This conclusion may however quickly change with more data. Stronger constraints on the Higgs invisible width will further constrain the resonance region and the combination of these constraints with future direct detection experiments may soon rule out the vector model.
Our study clearly demonstrates that, in the absence of positive signals, models of weakly-interacting DM particles will only remain viable if direct detection constraints can be systematically suppressed. This makes it increasingly interesting to study DM models with momentum-dependent scattering cross-sections. A systematic study of such theories will be left for future work. Conversely, Higgs portal models provide a natural framework for interpreting signals in the next generation of direct and indirect detection experiments. An advanced framework for such a reinterpretation using Fisher information will be implemented in future versions of GAM-BIT.
topes. Additional files can be provided to implement additional isotopes, and existing files can be replaced to study form factor uncertainties.
Of course, it is still possible to specify the WIMP coupling structure in the traditional way, e.g. by providing the effective couplings for spin-independent (SI) and spin-dependent (SD) interactions with protons and neutrons. In this case, DDCalc 2.0.0 will by default use the conventional form factors (i.e. the Helm form factor for SI interactions and the form factors from Ref. [149] for SD interactions, which can be found in DDCalc/data/SDFF/). In order to use the form factors from Ref. [147] also for standard interactions, one can set the global option PreferNewFF contained in DDCalc/src/DDConstants.f to true.
Let us finally emphasize that for general non-relativistic operators, the differential event rate depends not only on the conventional velocity integral f (v)/v d 3 v but also on the second velocity integral v f (v) d 3 v. As before, these velocity integrals are by default evaluated using the Standard Halo Model (SHM) with parameters that can be set externally. It is however also possible to provide tabulated velocity integrals in order to study velocity distributions that differ from the SHM. An illustration of this feature is provided in the example files DDCalc/examples/DDCalc_exclusionC. cpp and DDCalc/examples/DDCalc_exclusionF. f90, which demonstrates how to calculate an exclusion limit for a given WIMP model and a given velocity distribution.

A.2: Extended detector interface
The need to implement increasingly complex direct detection experiments has led to substantial extensions of how experiments can be defined in DDCalc 2.0.0. The details of this new interface are described in DDCalc/src/ DDDetectors.f, but we review the most important new features here.
First of all, it is now possible to define a number of different signal regions for each experiment and to specify the number of observed events and expected background events for each signal region. The simplest application is the implementation of a binned analysis, but it is also possible to define more complex signal regions, provided they can be characterised by a simple acceptance function (E R ), which quantifies the probability that a nuclear recoil with physical recoil energy E R will lead to a signal within the signal region. DDCalc 2.0.0 will then determine the expected signal in each signal region and calculate the binned Poisson likelihood. If the expected background in a signal region is set to zero, DDCalc 2.0.0 will interpret this to mean that the background level is unknown. In this case, the bin will only contribute to the total likelihood if the predicted number of signal events exceeds the number of observed events. The example files DDCalc/examples/DDCalc_ exampleC.cpp and DDCalc/examples/DDCalc_ exampleF.f90 illustrate how the predicted number of events in each signal region, as well as the resulting likelihoods, can be evaluated for specific parameter points.
Alternatively, one can also analyse experiments with unknown backgrounds using the optimum interval method by specifying the bins in such a way that their boundaries correspond to the energies of the observed events. Note that this means that it is typically not possible to use the binned Poisson method and the optimum interval method for the same choice of binning. A user wishing to compare these two analysis strategies should therefore implement them as separate experiments.
A related new feature is that it is now possible in DDCalc 2.0.0 to specify separate efficiency functions for each element (or indeed each isotope) in the target material. This is necessary for example if the efficiency of analysis cuts depends on the type of recoiling nucleus (as in CRESST) or if the low-energy threshold differs for different elements (as in PICO). For experiments with several different elements and several different signal regions, the number of efficiency functions that need to be specified can potentially be quite large. The preferred way to specify efficiency functions in DDCalc 2.0.0 is to provide external files with tabulated values, which by default are stored in DDCalc/data/. An illustration of this new structure can be found in the definition of the CRESST-II experiment (see below).
It is important to emphasize that the grid used to define the efficiency functions is also used to evaluate the other contributions to the differential event rate (i.e., form factors and velocity integrals). The number of grid points used in the definition of the efficiency functions directly influences the computation time and the precision of the result. In particular, it is essential to also provide a sufficiently large number of grid points in energy ranges where the efficiency is approximately constant. The function Retabula-teEfficiency in DDCalc/src/DDDetectors.f can be used to generate a fine efficiency grid from a coarse one, using linear interpolation between the provided values.

A3: New experiments
DDCalc 2.0.0 ships with a broad range of new experimental analyses. In particular, there are now a number of lowthreshold experiments, so that DDCalc 2.0.0 can now also be used to reliably calculate constraints on light DM. Moreover, we have implemented a number of planned experiments, which can be used to derive projected sensitivities.

CRESST-II
The CRESST-II results [91] are based on 52.2 kg days using the Lise detector module. Our implementation follows Refs. [150,151], i.e., we assume an energy resolution of σ E = 62 eV and take the cut survival probabil-ities from Ref. [150]. To avoid unnecessarily fine binning in energy ranges where the expected signal rates are small, we divide the energy range between 0.3 and 5.0 keV into 10 bins of increasing size. In the absence of a background model, we treat all observed events as potential signal events.
CDMSlite The analysis of CDMSLite is based on an exposure of 70.14 kg days [90,152]. The energy-dependent signal efficiency is taken from Ref. [90], which also describes the procedure for converting nuclear recoil energies into electron equivalent energies (eVee). We follow the same approach as in Ref.
[151] to determine the detector resolution, divide the energy range from 60 to 500 eVee into 10 bins of increasing size and assume no background model.

DarkSide-50
We implement the results from a search for heavy DM particles in the DarkSide-50 detector based on a total exposure of 19.6×10 3 kg days [93], taking the energydependent acceptance function from Ref. [93].
PandaX-II Since the most recent data taking period of the PandaX-II experiment (Run 10) has substantially lower background levels than previously analysed data sets [87,88], we implement it as an independent experiment (called Pan-daX_2017) rather than simply combining all runs. We use the same detector efficiency for the new data set as for our previous implementation of PandaX-II (see Ref. [101]) and assume a background expectation of 1.55 events. 15 It is then straight-forward to perform a combination of the different data sets by multiplying the individual likelihood functions.
XENON1T We use the same implementation of XENON1T [89] as described in detail in Ref. [79]. To reduce background levels, we focus on the central detector region with a mass of 0.65 t, and consider only events between the median of the nuclear recoil band and the lower 2σ quantile. We furthermore divide this signal region into two energy bins, which correspond to S1 ∈ [3 PE, 35 PE] and S1 ∈ [35 PE, 70 PE]. We estimate the expected backgrounds in the two bins to be 0.46 and 0.34 events, respectively, compared to 0 and 2 observed events.
LZ Our implementation of the LZ experiment [153] follows Ref. [154]. In particular, we assume an exposure of 5.6 · 10 6 kg days with a resolution of σ E /E R = 0.065 + 0.24 (1 keV/E R ) 1/2 and an acceptance of 50% for nuclear recoils. We consider 6 evenly-spaced bins in the range from 6 to 30 keV and assume a background of 0.394 events per bin. employ a C 3 F 8 target with 250 L fiducial volume. Six livemonths of data will be taken with a low threshold of 3.2 keV, which we implement using the same acceptance function as for PICO-2L [156], while 12 live-months will be taken with a threshold of 10 keV. We treat the two thresholds as two separate bins, in which case the expected backgrounds are 3 and 0.85 events, respectively. DARWIN The DARWIN experiment aims for a total exposure of 7.3 · 10 7 kg days with 30% acceptance for nuclear recoils and 99.98% rejection of electron recoils [157]. We assume an energy resolution of σ E /E R = 0.05 + (0.05 keV/E R ) 1/2 [158] and consider 5 equally-spaced bins between 5 and 20 keV. The dominant background is due to coherent neutrino-nucleus scattering, which we estimate from Fig. 3 in Ref. [158].

DarkSide-20k
We assume a total exposure of 3.65 · 10 7 kg days and estimate the energy resolution to be σ E /E R = 0.05 + (2 keV/E R ) [159]. To model the detector threshold, we implement the acceptance function for the f 200 -cut from Fig. 92 in Ref. [159]. We divide the energy range between 30 and 80 keV into 10 equally-spaced bins, and assume a background of 0.04 events per bin from instrumental background, as well as a total of 1.6 events (with non-trivial energy dependence) from coherent neutrino scattering, which we obtain by rescaling the results from Ref. [160].
Note that the number of observed events in each bin must be an integer in DDCalc, so it is typically not possible to set the observed number of events equal to the expected number of events in order to calculate the expected sensitivity of a future experiment. By default, the observed number of events is set to the integer closest to the background expectation, but this introduces a bias for example if there is a large number of bins with less than 0.5 expected background events. To accurately calculate expected sensitivities, one should simulate Poisson fluctuations in each bin, calculate the corresponding exclusion limits, and then construct the median exclusion. For an alternative approach, using Fisher information, we refer to Ref. [161].
Lastly, in Fig. 13 we show a comparison of the upper bounds on the spin-independent scattering cross-section determined using DDCalc with the official limits obtained by the respective collaborations. In all cases we find good agreement, validating our implementions of the experimental likelihoods in DDCalc. Also for the planned experiments described earlier we have confirmed that our sensitivity estimates are in sufficient agreement with the expectations published by the collaborations [153,155,157,159]. For fermion final states, the annihilation cross-section is given by where C f is a colour factor. For leptons, C f = 1, whereas for quarks, it includes an important 1-loop vertex correction given by [162] C f = 3 1 + For the hh final state, additional contributions appear from the 4-point contact interaction as well as DM exchange in tand u-channels. The annihilation cross-section for V V → hh is where the dimensionless quantities β = (1 − 2x h )/(v h v V ) and x = h m h /s, and v h and v V are the lab-frame velocities of the Higgs and vector DM, respectively. Similarly, the annihilation cross-section for χχ → hh (and equivalently for χ ↔ ψ) is given by