Global fits of simplified models for dark matter with GAMBIT

Global fits explore different parameter regions of a given model and apply constraints obtained at many energy scales. This makes it challenging to perform global fits of simplified models, which may not be valid at high energies. In this study, we derive a unitarity bound for a simplified vector dark matter model with an s-channel vector mediator and apply it to global fits of this model with GAMBIT in order to correctly interpret missing energy searches at the LHC. Two parameter space regions emerge as consistent with all experimental constraints, corresponding to different annihilation modes of the dark matter. We show that although these models are subject to strong validity constraints, they are currently most strongly constrained by measurements less sensitive to the high-energy behaviour of the theory. Understanding when these models cannot be consistently studied will become increasingly relevant as they are applied to LHC Run 3 data.


Introduction
As successful a theory as the Standard Model (SM) has been, there are many reasons for expecting it to exist within an even more descriptive particle theory.One of these reasons for beyond-Standard Model (BSM) physics is a number of astrophysical and cosmological observations that may require additional unseen matter [1][2][3].The WIMP hypothesis postulates that this matter consists of a Weakly-Interacting Massive Particle, and is a popular theory as it may explain the observed cosmological relic abundance of dark matter (DM) [4] and be strongly constrained by near-future experiments [5].
WIMP candidates are present in many UVcomplete theories including supersymmetric and extradimensional models.Rather than focus on these UVcomplete theories, this study will instead focus on a simplified model.These are a class of effective theories where the particle that mediates interactions between DM and SM particles is explicitly included.In the limit of large mediator masses, the traditional DM effective theory is recovered.These models have been reviewed in detail in many works, including Refs.[5][6][7][8][9][10][11][12].They have become the preferred method for modelling the simultaneous impact of low and high energy probes [13][14][15].Studies of these models are often grouped to include multiple simplified models with different mediator and DM spins.This work will instead focus on a single model, in which a vector DM candidate interacts with a vector mediator in the s-channel.Details of this model are discussed in section 2. For global fits of models with scalar or fermion DM candidates, we refer the reader to the previous work in this series [16].
Models containing new vector particles can come with additional theoretical challenges in the high energy limit of the theory, arising from the requirement of unitarity of the scattering matrix.Unitarity violation is a sign that the theory must be extended for it to be theoretically consistent; for example, unitarity violation in SM gauge boson interactions gave one of the early theoretical limits on the mass of the Higgs boson [17].Likewise, unitarity arguments have been used to place an upper bound on the mass of DM particles that obtain their relic abundance through thermal freeze-out [18].
Vector DM simplified models have been studied in detail for both high and low energy experiments.For direct detection constraints, it has been shown that additional non-relativistic effective operators may arise in these models [19,20], and that the use of polarized targets may distinguish between fermion and vector DM candidates [21].Assuming a detection of signal events at the XENONnT experiment, prospects for finding these models during Run III of the Large Hadron Collider (LHC) in dijet searches [22] and mono-jet searches [23] have been studied along with relic density limits [24].
In this work, we derive a unitarity bound from the self-scattering of vector DM and show the similarity in constraint between this and the requirement of a physical decay width of the mediator.We follow this with a global fit of this model using GAMBIT v2. 4, including the decay width and unitarity requirements in our calculations.This paper is structured as follows.In Section 2, we describe the simplified model that we study, and the reasons behind the choice of couplings.In section 3, we derive a unitarity bound on this model.Section 4 describes the set of experimental constraints we use to perform a global fit of this model and section 5 provides our results.Finally, Section 6 briefly discusses the potential to observe these particles at near-future experiments and presents our conclusions.The samples from our scans, the corresponding GAMBIT plotting scripts and a detailed unitarity bound proof can be downloaded from Zenodo [25].

Model
The general form of the Lagrangian for a simplified model of vector DM X µ coupled to quarks via a mediator V µ with vector and axial-vector couplings is [23] where X µν is the field strength tensor for the vector DM, and F ′ µν for the mediator.To reduce the complexity of this simplified model and the dimensionality of the corresponding parameter space, we make a number of simplifying assumptions.First, we neglect any four-field interactions, which are expected to be irrelevant for phenomenology, and therefore set the couplings λ DM , λ M , b 3 and b 4 to zero.Furthermore, we assume that the simplified model conserves CP symmetry, which requires the real components of b 6 and b 7 in eq. ( 1) to vanish.Finally, to preserve the SM gauge structure, we concentrate on vector-like couplings of the mediator to SM quarks and set h 4 = 0.
With these restrictions, one finds that the imaginary components of b 6 and b 7 only give rise to interactions that vanish in the limit of zero momentum transfer, leading to strongly suppressed constraints from direct detection experiments.Including these couplings in our global fits would therefore lead to rather trivial results, while at the same time requiring significant additional work in order to correctly treat the non-relativistic effective operators O 19 and O 20 introduced in Ref. [20] and the interference between different operators in the simulation of LHC events.We therefore neglect these couplings in the present work and focus on the two interaction terms proportional to h 3 and b 5 .
Therefore, the Lagrangian of the model we adopt is where we choose to label the quark coupling as g q and the DM coupling as g DM to agree with our previous work [16].Both couplings can be taken as purely real since any imaginary phase can be absorbed into a redefinition of the fields.Perturbative unitarity breaks down in large regions of the parameter space of this model due to the poor high energy behaviour of the longitudinal polarized modes of the vector DM.Following the same approach as Ref. [26], here we derive an approximate unitarity bound for this model in terms of the Mandelstam variable s, from scattering of vector DM Section 3 derives this relation, and section 4.4 describes how unitarity was imposed on simulated collider events in our global scan.In Appendix A, we present the equivalent bound if the b 6 and b 7 couplings of eq. ( 1) are included alongside the b 5 coupling.
The onshell decay width of the mediator to a pair of DM particles, V → XX, is and the width to a given pair i of SM quarks, V → q i q i , is The total width of the mediator should not exceed the mediator mass, or else the perturbative description of DM interactions via mediator exchange is expected to break down.

Forming Unitarity constraints from partial waves
Unitarity bounds are formed from partial wave analysis of the scattering of vector DM particles.For examples on the use of this method, see e.g.Refs.[17,26,27].From the requirement of partial wave unitarity, the scattering amplitude must obey the bounds and Here M J ii is the full scattering matrix element between 2particle states i where the initial and final state particles are the same (hence the repeated index i), for the Jth partial wave.Tree-level amplitudes are generally used to form these bounds, assuming that the higher orders do not provide significant corrections to the amplitude.In this way, the resulting bound may be interpreted as a "perturbative unitarity" bound.In the case of zero initial and final total spin, Here P J (x) is the Legendre polynomial of order J, θ is the scattering angle and s is the square of the centreof-mass energy.An additional factor of 1/ √ 2 must be applied to the right hand side for each initial or final state with identical particles.The term β ii is a kinematic factor, which for a final state of equal mass DM particles becomes In the high-energy limit (s → ∞), β ii approaches 1.As the zeroth order usually dominates, it is often sufficient to study In the following derivation, we consider the selfscattering of DM, rather than DM with its antiparticle.The particle-antiparticle scattering via s-channel mediator exchange will also face poor behaviour at high energies, however this will be effectively covered anyway by our additional requirement that the perturbative description of the off-shell decay width of the mediator (including to DM particle-antiparticle pairs) does not break down.The tree level amplitude of DM self-scattering has contributions from t and u channel processes (see Figure 1), which can be derived separately, and summed together.This is most easily understood in the centre of mass frame, where for incoming particles (with momenta p (1) and p (2) ) and outgoing particles (with momenta p (3) and p (4) ), p (1) = E, 0, 0, P p (2) = E, 0, 0, −P p (3) = E, P sin θ, 0, P cos θ p (4) = E, −P sin θ, 0, −P cos θ .(11) Here E = Ecm 2 is the incoming particle energy and P is the magnitude of the incoming momentum of each particle.The longitudinal polarisations will most strongly violate unitarity, and so it is sufficient to solely form a bound from evaluating the amplitude for incoming longitudinally polarised DM particles.In the centre of mass frame, these are The amplitude for t-channel DM-DM scattering at treelevel is where Evaluating this amplitude in the centre of mass frame gives Similarly, the scattering amplitude for u-channel DM DM scattering at tree-level is

Unitarity Bound
The total amplitude of the scattering process is Performing the integral in eq. ( 10) and substituting into eq.( 7) gives the bound on parameters to satisfy unitarity Since unitarity is increasingly violated as the collision energy increases, the limit s ≫ m 2 DM is often taken in the literature.If this limit is taken, this bound simplifies to The validity of this limit breaks down for small DM masses and large couplings.In these cases, the complete bound eq. ( 18) should be used.Even though the unitarity requirement above has been derived for the case of DM self-scattering, the resulting bound can be interpreted more generally as the energy scale where the interactions between DM particles and the vector mediator become unphysical.We will therefore apply the unitarity bound from eq. ( 18) to any process in which a pair of DM particles is produced, with √ s being replaced by the invariant mass of the DM pair m inv .In particular, this requirement will be implemented in our simulation of LHC monojet events (see section 4), where we will discard any event that violates the unitarity bound.In other words, we apply LHC constraints only on those regions of phase space where the simplified model predictions can be trusted, and set conservative bounds otherwise.
It is worth noting that for m DM < m M /2, we expect mono-jet production to proceed dominantly via an onshell mediator, such that m inv ≈ m M .Hence, for virtually all events will be removed by the unitarity requirement such that the LHC mono-jet bounds are effectively absent.However, parameter points in this region typically also violate the requirement on the decay width from eq. ( 21), such that they would be excluded from the analysis anyway.

Physical Decay Widths
Alongside unitarity violation, another indication that the model breaks down is that the decay width of the mediator becomes large, indicating the inapplicability of perturbation theory to the underlying scattering process.When the mediator is on-shell, this can be interpreted as a bound on the decay width We reject all points in parameter space that do not satisfy this bound.In the following we require that an analogous inequality also holds for the off-shell decay width when replacing m M by √ s: In the high energy limit, the bound on the off-shell decay width results in the requirement This differs from eq. ( 18) by a factor of 2 (the unitarity bound being the stricter of the two).When assuming high collision energies, it is therefore clear to see that the unitarity bound and off-shell decay width bound are practically interchangeable.Figure 2 shows a comparison between the unitarity constraint with and without The requirement of a physical off-shell decay width (red) excludes a smaller region than the requirement of unitarity (green), but follows a similar trend.Taking the high-energy limit of the unitarity bound is a consistently stricter cut on the parameter space (blue).taking the high-energy limit, for a representative choice of parameters, along with the exclusion from requiring that the off-shell decay width is physical.The similarity between the unitarity and decay width conditions would suggest that for the choice of parameters shown, very little difference would be observed if the two were interchanged.

Constraints
Interactions between DM and SM quarks are constrained by many different measurements of astrophysical, cosmological and particle physics processes.
We use likelihoods, implemented in GAMBIT 2.4, for DM direct and indirect detection experiments, collider searches at the ATLAS and CMS experiments, and the measurement of the DM relic abundance.We generate the necessary model-specific GAMBIT module functions (including those used to store spectrum and decay information [28]) using the GAMBIT Universal Model Machine (GUM) [29].This includes interfaces to backend codes that contain physics calculations for each DM observable.We apply the perturbative unitarity and physical off-shell decay width constraints described in section 3. limits are conservative; this is detailed in section 4.4.
We reject parameter points that fail the requirement of a physical on-shell decay width of the mediator, before calculating their likelihood contributions.
Table 1 provides a summary of each likelihood that we include that is sensitive to BSM physics.For each likelihood, we provide either: ln L bg , the value that the likelihood takes purely from the SM, or ln L max , the bestcase likelihood that can be achieved when parameters exactly match their centrally measured values.
For a detailed description of the implementation of each likelihood in GAMBIT, we refer the reader to the previous work in this series [16].We provide brief summaries of each likelihood in the following subsections.

Relic Density
We use GUM to generate the CalcHEP v3.6.27[55,56] model files that are supplied to micrOMEGAs v3.6.9.2 [57].The relic density of DM is obtained with the DarkBit interface which uses micrOMEGAs to solve the Boltzmann equation for the number density of DM particles in thermal equilibrium, assuming a standard cosmological history.To form a likelihood from the relic abundance, we compare the calculated density to the Planck 2018 measurement of Ω DM,obs h 2 = 0.120 ± 0.001 [54] with a 1% theoretical error added in quadrature with the quoted Planck uncertainty.
We study both cases where the DM candidate is a subcomponent of the observed relic abundance and where it fully saturates the abundance.When requiring that it saturates the relic abundance, we use the Planck measurement to define a Gaussian likelihood based on the predicted WIMP abundance.When allowing it to form a subcomponent, we modify this likelihood to be flat for predicted densities below the measurement; details can be found in Ref. [58].

Direct Detection
The parameters of a simplified DM model can be translated to the coefficients of the relevant operators in a non-relativistic EFT for WIMP-nucleon scattering, c N i (q 2 ).The single relevant operator and its coefficient for the vector DM simplified model in this study is [23].
which was supplied to DDCalc v2.2.0 [59,60], to compute the differential cross-section and target element of interest.We do not include the effect of operator mixing from running as it has been shown to have little effect for pure vector couplings of the mediator to quarks [13].

Indirect Detection
The model we study has two primary DM annihilation channels, annihilation to mediators and to quarks.Annihilation to a pair of mediators occurs as an s-wave process, and will be the primary annihilation channel when kinematically allowed (m DM > m M ).When this channel is closed, the annihilation will occur to a pair of quarks, through the suppressed p-wave channel.We do not include p-wave contributions to the gamma-ray flux as they should not be large enough to impact searches toward dwarf spheroidals for the model we consider.
We compute the annihilation cross-section with CalcHEP, using the GUM interface to generate the required CalcHEP model files.We use the combined analysis of 15 dwarf spheroidal galaxies, Pass-8, performed by the Fermi-LAT Collaboration over 6 years of data taking [53], using gamLike v1.0.1 to compute the likelihood through its interface to DarkBit.DM annihilations at the centre of our galaxy are an alternative to dwarf spheroidal measurements.Since Fermi-LAT Galactic Centre limits are not as robust as limits from dwarf spheroidals, we do not include them in this study.We do however briefly comment on the future impact of CTA observations on the parameter space of this model in section 6.

Monojet searches at the LHC
One of the primary channels via which to search for the model at colliders is the creation of a pair of final state WIMPs in association with a jet created by initial state radiation.This gives a signature of a single jet plus missing transverse energy ( / E T ).We include the most current monojet searches from CMS and ATLAS searches with 137 fb −1 [52] and 139 fb −1 [51] of Run II integrated luminosity respectively.
To calculate the total production cross-section σ and the product of the efficiency and acceptance for passing the analysis kinematic selections ϵA, we perform simulation of Monte Carlo events with MadGraph_ aMC@NLO [61] (v3.1.1),interfaced to Pythia v8.3 [62] for parton showering and hadronization.To form the quantity ϵA we pass these events through MadAnalysis 5 [63] and implement the ATLAS and CMS monojet analyses.Rather than perform this calculation for each parameter sample, we precompute a grid of cross sections (σ) and ϵA factors in advance, and interpolate them at runtime using ColliderBit [64].
An additional analysis cut is added to our implementations of the ATLAS and CMS kinematic selections, to remove any events which would violate the unitarity bound presented in section 3, replacing √ s with the invariant mass of the DM pair.When this cut becomes strong enough, there is a significant drop in the predicted acceptance of the analysis, and we can no longer make any sensible predictions regarding collider constraints.If no simulated events pass the unitarity cut, we expect the parameter point to be unobservable at the LHC and simply assign the background-only likelihood.
The interpolation grid we use is as follows: - The grids for the mediator mass and couplings were chosen to be approximately equally spaced in log-space.The ratio of DM and mediator masses is more effective than the DM mass as a grid variable as it allows us to choose a grid with a higher density of points across the resonance region, where we expect rapid changes in predictions.Below the DM mass/mediator mass ratio of 0.01, we assume that we can safely extrapolate to small DM masses as the predicted signal should not vary significantly.After removing any points with DM masses above the limits of our scan, this gives a total number of 6370 grid points.

Searches for dijet resonances
The presence of a mediating particle in the model may generate dijet events at colliders, with an invariant mass of approximately the mediator mass.Dijet resonance searches provide robust constraints on DM simplified models, where the extremely high multijet background must be removed with clever kinematic analysis cuts.
The cross-section for the production of a dijet resonance can be approximated as the product of the crosssection of mediator production and the branching ratio of the mediator into quarks, assuming that the narrow width approximation holds.When the ratio of the mediator decay width to mass is high, this approximation breaks down, and our treatment of dijet searches would become dubious.We briefly investigate the dependence of the model exclusion on this assumption in section 5.
We implement dijet limits provided by ATLAS and CMS [42-50] by scaling of the published limits of the mediator-quark coupling by the branching ratio into quarks, following the same approach as Refs.[16,65].These published limits are interpolated in m M for each parameter point, and the likelihood is formed from the most constraining search for a given mediator mass.The combined coupling upper limits are provided in Figure 1 of Ref. [16].
In the absence of tree-level couplings of the mediator to leptons, couplings at loop level may still be generated through kinetic mixing, and the model may be observable at dilepton searches.Despite the tight constraints on dilepton signatures for vector mediated simplified models, the loop suppression of these lepton couplings will prevent dilepton constraints on the quark coupling being any stronger than dijet limits.For this reason we do not include dilepton constraints in this study.For a discussion on the lepton couplings generated through kinetic mixing, we refer the reader to Refs.[26,66].

Nuisance Parameter Likelihoods
Along with the model parameters in the model we study, we also include a set of nuisance parameters which are used in each of our astrophysical likelihoods.A complete list of these parameters is given in Table 2.
We treat the local DM density ρ 0 following the standard procedure in DarkBit, where ρ 0 is assumed  to be log-normally distributed, centred around ρ 0 = 0.40 GeV cm −3 and with an error σ ρ0 = 0.15 GeV cm −3 .The scan range of ρ 0 is asymmetric to reflect this distribution.3σ ranges for all other nuisance parameters are provided in Table 2.
We treat the Milky Way halo in the same way as in several of our previous DM studies [16,60,67], where the DM velocity is assumed to follow a Maxwell-Boltzmann distribution.The peak velocity and Galactic escape velocity uncertainties are described by Gaussian likelihoods with v peak = 240 ± 8 km s −1 [68] and v esc = 528 ± 25 km s −1 (based on Gaia data [69]), respectively., where the ideal likelihood is the combination of background-only and maximum possible likelihoods detailed in Table 1.

Results
We have performed a comprehensive scan of the model parameter space using the differential evolution sampler Diver v1.0.4 [70] with a convergence threshold of 10 −6 and a population of 20 000, with an additional scan for DM masses below 2 TeV to improve sampling.We carried out two separate scans for the case where the observed DM relic density is taken as an upper limit or as a two-sided measurement.Unlike the previous study in this series [16], scans with a capped LHC likelihood were not performed, as any small preferences over the background-only hypothesis in mono-jet searches were not found to occur within the surviving parameter space of the scan.A scan with a capped LHC likelihood would therefore produce results that were indistinguishable from its uncapped equivalent.
Table 2 provides the full list of parameters and scan ranges.We adopt the same choice of scan ranges and sampling distributions of the masses and couplings as those in Ref. [16].Very small couplings are avoided in order to focus on regions where unitarity violation may be relevant.The coupling upper bounds are of order unity in order to keep the decay width of the mediator from becoming excessively large.The range of masses was chosen to focus on regions where it was expected that both collider searches and direct and indirect searches may be complementary.The parameter points that give the best likelihoods are given in Table 3.
The profile likelihood from combined constraints on the complex vector DM model is shown in Figures 3  and 5.The model prefers parameter regions where DM annihilation is efficient, and there are two regions corresponding to the two DM annihilation channels.Around the diagonal m M ≈ 2m DM , the annihilation occurs close to a resonance into a pair of quarks.For regions where m DM > m M , the annihilation occurs as a t-channel process into a pair of mediator particles.Below approximately 500 GeV, the annihilation may not be great enough to prevent exclusion from direct detection constraints without leaving the limits of the scanned parameter ranges.This shape is highly similar to those presented for a scalar DM candidate in Ref. [16].This is because the strongest limits come from the direct detection experi-ments, which are dependent on the effective operators that are relevant, and this model shares the same relevant operator as the scalar DM model.The model survives for a greater proportion of the parameter space than the scalar DM model, despite the additional inclusion of PandaX-4T direct detection data in this work.The small variation in the profile likelihood around 2 TeV is a sampling artifact, and does not reflect any physical change in predictions.
In Figure 4, we show how the profile likelihood changes if the scan range were extended to masses up to 100 TeV.The resonance region closes off around 30 -40 TeV as DM becomes overabundant unless the couplings become non-perturbative.The non-resonant region continues on with a largely flat likelihood.For DM masses well beyond 100 TeV, thermally produced DM will violate generic unitarity bounds [18].
Requiring that the DM relic abundance is saturated shrinks the surviving region to mediator masses above 1 TeV for the off-resonance region.For lower mediator masses, the non-relativistic effective coupling to nucleons is stronger and therefore expected signal at direct detection experiments is greater.Figure 6 (left) shows that at low mediator masses, the likelihood is higher in parameter regions where the model strongly underproduces DM to avoid tension with these experiments.As the strength of the direct detection constraints increases toward lower DM mass, the surviving parameter region also has a lower bound on the DM mass that may be seen in Figure 6 (right).The surviving region along the resonance does not depend strongly on whether the abundance likelihood is taken as a one-sided upper limit or a two-sided measurement.Measurements of dwarf spheroidal galaxies do not appear to have any strong influence on the profile likelihoods.
We find that, in the surviving parameter regions, the decay width of the mediator is dominated by the partial width to quarks.Limits from dijet searches prevent mediator-quark couplings g q above roughly 0.1 for most of the parameter space.This preference toward lower g q reduces the effect of high decay widths, as the partial width to quarks is proportional to g 2 q .Fig 7 shows that within 2σ of the best-fit point, the decay width of the mediator does not exceed 0.02m M , safely satisfying the narrow width requirement.The effect of monojet searches cannot be seen directly on the results of the profile likelihood.For any model parameters where monojet searches would have sensitivity, these are strongly excluded by relic abundance limits and direct detection searches.The combined global fit therefore does not appear to be strongly affected by unitarity considerations.This conclusion might however change when considering a more general parameter space including also the couplings b 6 and b 7 .
The best fit for each scan lies along the resonance, at the upper limits of the masses, and toward the lower limits of the quark coupling.In these regions, the relic abundance, and the strength of the direct detection signals are minimised.When the DM candidate is allowed to be a subcomponent of the observed DM density, this best fit point approximately matches the background likelihood as the signals at any given DM experiment are almost entirely negligible.We compute an approximate p-value of the best-fit likelihood conditioned on the 'ideal' scenario (sum of background-only and max entries in Table 1) for 1-2 effective degrees of freedom.Further explanation of the construction of this particular p-value can be found in Ref. [71].Neither case (saturated or subdominant DM) is disfavoured, returning p-values of 0.3 and above.
We limited the couplings to be no lower than 0.01, in order to target parameter regions where unitarity violation was most likely to cause issues without introducing large hierarchies between couplings.If the scan range was expanded to include smaller g q , it can be seen from Figure 5 how the size of the surviving parameter space should increase.Expanding the lower limit on g DM will only expand the surviving space if g q is also expanded.For the parameters scanned over in this work, the model is excluded for lower g DM , as there cannot be sufficient annihilation of the thermal DM abundance.

Discussion
In this work, we have derived a unitarity bound for a simplified model with a vector DM candidate that interacts with SM quarks via an s-channel vector mediator.We showed that this unitarity bound is highly similar to the bound on the model parameters one would require from the behaviour of the off-shell decay width, which is another challenge that plagues these theories.Applying this bound to simulated collider events, we performed a global scan of this model with GAMBIT.We found that in all of the simulated parameter regions where the unitarity of the model may come into question or the decay of the mediator becomes unphysical, the model is excluded by experiments that are less sensitive to the high energy behaviour of the theory.Since the model exclusion most strongly comes from direct detection experiments and relic abundance limits, the surviving parameter space is split in two by the DM annihilation channels.The overall result is a series of limits that are highly similar to, but slightly weaker than, those found for corresponding scalar and fermionic DM models in the previous study in this series [16].
In the coming years, many experiments are expected to take data that may be used to constrain the model that we consider.In Figure 8 we show the predicted number of signal counts at the next-generation liquid Xenon direct detection experiment, DARWIN [72].Within the surviving parameter space of the model, up to several hundred recoil events may be observed.Depending on how effectively the background can be rejected, a large portion of the surviving parameter space in these scans may be ruled out in the absence of any signal measurements.
We also checked the extent to which future observations by the Cherenkov Telescope Array (CTA) would constrain the model, using the same methods as in Ref. [16].None of the currently viable parameter space will be probed by CTA, with the parameter space along the resonance region far out of reach because the annihilations occur through the p-wave suppressed channel to quarks.
Finally, we note that further constraints can be expected from Run 3 of the LHC and the subsequent highluminosity phase, as well as future colliders.In order to correctly interpret these constraints, it will become increasingly important to understand the high-energy behaviour of simplified models.1) are allowed to be nonzero, the unitarity bound becomes  The term from the real component of the b 7 coupling is independent of s.

Fig. 2 :
Fig. 2: Comparison between unitarity violation and unphysical decay widths for a demonstrative choice of parameters (s = 10 8 GeV 2 , m M = √ s, gq = 0, and varying m DM and g DM ).The requirement of a physical off-shell decay width (red) excludes a smaller region than the requirement of unitarity (green), but follows a similar trend.Taking the high-energy limit of the unitarity bound is a consistently stricter cut on the parameter space (blue).

★ΩFig. 3 :Fig. 4 : 8 ]
Fig.3: Profile likelihood, profiled over couplings.The measured DM relic abundance is taken as an upper limit (left) or to be composed entirely of the vector DM candidate (right).1σ and 2σ contours are shown in white, with the star representing the best-fit point.

Fig. 6 :
Fig.6: DM relic abundance for the surviving parameter space, when taking the relic abundance measurement as an upper limit.We show the abundance both against mediator mass (left) and against DM mass (right).1σ and 2σ contours are shown in white.

Fig. 7 :
Fig. 7: Profile likelihood, as a function of the mediator width to mass ratio, profiled over all model parameters.1σ and 2σ confidence limits are shown in black, with the red star representing the best-fit point.

ΩFig. 8 :
Fig. 8: Predicted number of signal events in the DARWIN experiment, coloured by the mediator mass.1σ and 2σ profile likelihood contours are shown in white.
in part performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk).The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1.DiRAC is part of the National e-Infrastructure.PS acknowledges funding support from the Australian Research Council under Future Fellowship FT190100814.TEG and FK were funded by the Deutsche Forschungsgemeinschaft (DFG) through the Emmy Noether Grant No. KA 4662/1-1 and grant 396021762 -TRR 257.MJW is supported by the ARC Centre of Excellence for Dark Matter Particle Physics (CE200100008).This article made use of pippi v2.2 [73].A: Unitarity Bound including b 6 and b 7 couplings If the b 6 and b 7 couplings from eq. (

M m 2 M−Re(b 7 ) 2 m 4 DM s 2 b 5
corresponds to the coupling g DM in the model we adopt.For the proof of the relation, we refer the reader to the supplementary Zenodo record for this study[25].The b 6 and b 7 couplings are split into their real and imaginary components, with the CP-violating couplings left in for completion.The imaginary component of the b 7 cancels in the formation of the bound.In the limit of high s, this simplifies to

Table 1 :
2 to the calculation of collider signals, to ensure that calculations are accurate and the resulting All likelihoods included in our fits.We give the SMonly (i.e.background-only) log-likelihood ln L bg for those that search for events above an SM background.For the rest, we give the highest achievable value of the log-likelihood ln L max , where the predicted value of the chosen observable or a nuisance parameter is exactly equal to its measured value.

Table 2 :
List of model and nuisance parameters and their corresponding scan ranges.

Table 3 :
Relic Density Best Fit m DM (GeV) Best Fit m M (GeV) Best Fit gq Best Fit g V Approximate best-fit points for each scan.∆ ln L values are defined as ln L − ln L ideal