Di-jet/$e^+e^-$+ MET to Probe $Z_2-$Odd Mediators to the Dark Sector

We explore a scenario where Dark Matter (DM) couples to the Standard Model mainly via a scalar mediator ${\cal S}$ that is odd under a $Z_2$ symmetry, leading to interesting collider signatures. In fact, if linear interactions with the mediator are absent the most important DM production mechanisms at colliders could lead to final states with missing transverse energy (MET) in association with at least two fermions, such as di-jet or di-electron signatures. The framework we consider is model-independent, in a sense that it is only based on symmetry and formulated in the (extended) DM Effective Field Theory (eDMEFT) approach. Moreover, it allows to address the smallness of first-generation fermion masses via suppressed $Z_2$ breaking effects. From a di-jet +MET analysis at the LHC, we find rather loose bounds on the effective ${\cal S}$-${\cal S}$-DM-DM interactions, unless the mediator couples very strongly to SM fermions, while a future $e^+ e^-$ collider, such as CLIC, could deliver tighter constraints on the corresponding model parameters, given the mediator is leptophilic. We finally highlight the parameter space that allows to produce the observed DM density, including constraints from direct-detection experiments.


I. INTRODUCTION AND SETUP
The origin of the dark matter (DM) observed in the universe is one of the biggest mysteries in modern physics. It is tackled by a multitude of experiments, which are currently running or in preparation and are probing very diverse energies. While experiments aiming for a direct detection (DD) of DM particles via nuclear recoil typically feature collision energies in the keV range, collider experiments, trying to directly produce DM particles, probe momentum transfers exceeding the TeV scale. Combining results from all such kinds of experiments in a single, consistent, yet general framework is important in order to resolve the nature of DM.
In [1], such a framework to describe and compare searches at different energies was proposed, based on effective field theory (EFT), however allowing for detectable collider cross sections without relying on the problematic high energy tail of distributions [2, 3] and reproducing the correct relic density while avoiding a (too) low cutoff. To this end, in the eDMeft approach [1], the field content was enlarged by a dynamical (pseudo-)scalar (and potentially light) mediator S to the dark sector, the latter being represented by a scalar or fermionic field χ. Since both the mediator and the DM are assumed to be singlets under the SM gauge group, they can in principle interact via renormalizable couplings, however fully consistent interactions of the mediator with SM fermions (or gauge bosons) require D = 5 operators due to gauge invariance -which are not incorporated in typical simplified DM models [4][5][6][7]. In the eDMeft such couplings are included properly in the EFT framework, which is then consistently truncated at the D = 5 order, leading to a well controllable number of new parameters and avoiding the need to stick to a specific UV completion. The inclusion of the most general set of (non-redundant) D = 5 operators, allows in particular to consider richer new physics (NP) sectors, than just consisting of a single dark state and one mediator.
In this paper we focus on the phenomenology of the D = 5 operator S 2χ L χ R , which can give rise to interesting di-jet phenomenology at colliders, as we will see below. If for example symmetries forbid the dimension four Sχ L χ R interaction, this coupling could in fact be the main portal to the Dark Sector, which could be missed in DD experiments, while mono-jet searches should be adjusted to take advantage of the peculiar di-fermion final state.

A. General Setup
We thus start from the effective Lagrangian of the SM field content, augmented with a fermion DM singlet χ and a real, CP even scalar mediator S, including operators up to D = 5, as presented in [1], with the additional assumption that the coefficient of the operator Sχ L χ R is constrained to be negligibly small. 1 For concreteness we will assume in the following that a symmetry forbids such D = 4 interactions with the DM, where the most simple choice is assuming S to be odd under a Z 2 parity, S Z 2 −−→ −S, under which we take all the SM fields to be even, with the exception of the right-handed first fermion generation, which is also odd.
Beyond entertaining a new portal to the dark sector which is testable at (future) particle colliders, yet in agreement with null-results in DD so far, this scenario can also motivate the smallness of first-generation fermion masses, which are now forbidden at the renormalizable level. 2 Eventually, many of the terms in this modified eDMeft vanish compared to the original setup [1], including those with an odd power of mediators, unless they feature the right-handed up or down quark (or the corresponding electron). On the other hand, as mentioned, the SM-like Yukawa couplings of the latter fermions vanish and the corresponding masses will thus only be generated via small Z 2 breaking effects equipped with cutoff suppression. The corresponding Lagrangian reads where Q L and L L are the left-handed SU(2) L quark and lepton doublets, resp., d R , u R , and R are the righthanded first-generation singlets, and H is the Higgs doublet. 3 The latter develops a vacuum expectation value (vev), | H | ≡ v/ √ 2 174 GeV, triggering electroweak symmetry breaking (EWSB). In unitary gauge, the Higgs field is expanded around the vev as H 1/ √ 2(0, v +h) T . Here, h is the physical Higgs boson with mass m h ≈ 125 GeV. Finally, L SM denotes the SM Lagrangian without the Yukawa couplings of the first generation, see Eq.
In contrast to the original setup, we assume the mediator to develop a small vev | S | ≡ v S ∼ O(1−10) MeV, which finally generates masses for the first fermion generation. Since the resulting mixing with the Higgs via the |H| 2 S 2 operator is suppressed, the latter will not be considered in the following. Finally, also the "usual" dark matter coupling Sχχ is generated by the spontaneous breaking of the Z 2 -symmetry, with coefficient ∼ 2y S χ v S /Λ, which is however highly suppressed and only plays a role in direct detection experiments, see below. The coefficient of the potential second D = 5 portal to the dark sector allowed by the symmetry, |H| 2χ L χ R , will on the other hand taken to be small from the start, as motivated to evade direct detection constraints (remember that v/v S ∼ O(10 4 )) and limits from invisible Higgs decays (for light dark matter) [8], playing therefore no role in the collider discussion.
Neglecting leptons for simplicity, which can be treated analogously, the resulting mass terms read (2) where q = u, d are three-vectors in flavor space and the Yukawa matrices reflect the Z 2 assignments. Without breaking of the latter symmetry via v S > 0, one quark family would remain massless, corresponding to a vanishing eigenvalue of Y H q . On the other hand, a small breaking of v S ∼ O(10) MeV is enough to generate appropriate m u ∼ m d ∼ 5 MeV with O(1) Yukawa couplings and Λ 1 TeV.
After performing a rotation to the mass basis with U d L = U u L V CKM , we obtain the couplings of the physical quarks to the Higgs boson and the scalar mediator where in particular the latter are crucial to test the S 2 χ 2 operator at colliders, relying on a coupling of the mediator to the SM.

B. Flavor Structure
To fully define the model, we need to fix a flavor structure, avoiding excessive flavor-changing neutral currents (FCNCs). The latter are generically generated since the fermion mass matrices M q receive contributions from different sources (see Eq. 2) and are in general not aligned with the individual scalar-fermion couplings ∼ Y H,S q , such thatŶ H,S q will not be diagonal. To this end, we first note that, in the interaction basis, the Yukawa matrices can be expressed in terms of the mass matrices as In the mass basis, they becomê where the unitary rotations of the left-handed fermion fields drop out since they share the same Z 2 charges and their couplings (with a fixed right-handed fermion) are thus aligned with the corresponding mass terms. This is not true for the right handed fermions, where the corresponding rotation matrices induce a misalignment and thus FCNCs. However, while it would not be possible to entertain U u L = U d L = 1, since then V CKM = 1, in conflict with observation, one can in fact choose the Yukawas matrices in Eq.  4 Although a more systematic analysis of FCNCs in such a scenario would be interesting, we will just stick to the latter choice for the rest of this article, ending up with only diagonal cou-plingsŶ This means that the second and third generation couple to the Higgs boson as in the SM while the first generation couples instead only to the DM mediator, with strength determined by the free parameter v S , which we will trade for y S u /Λ ≡ (Ŷ S u ) 11 /Λ in the following. While the latter should not be too tiny, since then a very large Z 2 -breaking vev v S will be required to reproduce the quark masses, as discussed, O(1) values of y S u v/Λ are in perfect agreement with a modest vev and a reasonable cutoff.
So far we did not include the lepton sector, however a similar setup is possible for the latter, leading straightforwardly toŶ Finally, expressing everything in terms of y S u , we obtain the relations for the couplings of the mediator to SM fermions, plugging in the values m u = 2.5 MeV, m d = 5 MeV, m e = 0.5 MeV. As mentioned, y S u /Λ can be chosen basically freely, however should not violate perturbativity of the EFT (and of the potential UV completion), which con- where we made use of the fact that the S−Yukawa scales like y S f ∼ g 2 UV .

C. Relevant Parameters
In the following, we will derive the prospects to constrain the Z 2 -symmetric bi-quadratic portal S 2χ L χ R and the S-Yukawa coupling from LHC and future (e + e − ) collider data, meeting constraints from DD and the observed relic density. A unique process where the new portal enters is fermion-pair-associated DM production, as induced by the Feynman diagrams given in Fig. 2, with the DM leading to a characteristic missing energy signature. Before moving there, we will however summarize the relevant physical parameters in the model at hand. These are • the bi-quadratic portal coupling y S χ /Λ • the S−Yukawa coupling y S u /Λ , where we neglected potential scalar mixing from λ HS . 5 While this defines the main model being studied in the following sections, there are also two interesting variants obtained by either assigning positive Z 2 parity to all leptons or to all quarks. This will lead to a leptophobic or hadrophobic mediator, respectively, with y S e = 0 and finite y S d = 2y S u or vice versa.

AT (HL)-LHC
To get a first idea on near-future constraints on the new DM portal, we derive bounds from current (and 5 In the following analysis, we will consider the mediator to be much heavier than its vev, which requires an additional contribution to the Lagrangian (1). While a cubic term needs a very large (non-perturbativ) coefficient, a straightforward possibility is to add another singlet S 2 , already envisaged in footnote 2, with a O(TeV 2 ) quadratic term and a mass mixing SS 2 with O(1 GeV 2 ) coefficient and/or a SS 3 2 portal with coefficient O(10 −6 ). We have checked that other effects of the new scalar can be effectively decoupled.
projected future) LHC runs employing the CheckMate implementations of existing ATLAS analyses. A unique signature to constrain y S χ is di-jet production in association with MET, see Fig. 2 with the electrons replaced by up or down quarks. Here, the new portal enters at the tree-level, while the main background is νν production in association with jets. Although a dedicated analysis on the particular di-jet topology could improve the sensitivity, we expect the existing mono-jet search [10] using 36.1 fb −1 of data and a SUSY motivated search for multiple jets plus missing energy [11] to deliver already relevant constraints. Thus, we refrain from setting up a custom analysis but rather focus on future leptonic colliders for that purpose, where in particular the large QCD backgrounds faced at the LHC are avoided and the limits are expected to be much stronger.
Regarding the mentioned LHC analyses, the latter one naively delivers stronger constraints, but here events are used that have energies above the envisaged cutoff Λ = O(1) TeV such that the validity is questionable [2, 3,12]. The scalar sum of the transverse momenta of the leading N jets and E miss T are required to be at least 1.6 TeV. Therefore a reasonable value for the cut-off is at least Λ 3 TeV. In addition all signal-regions are inclusive ones, which means that they include events with even much higher energies, such that the resulting constraints would only be valid for borderline large couplings y S u . Exclusive signal regions (EM), as provided in [10], allow for a better estimate of the event energy. For that reason we constrain ourselves to signal regions up to EM6 of [10], the latter containing events with E miss T = (600 − 700) GeV, to get robust constraints.
The actual bounds on the couplings and the prospects for the HL-LHC with a luminosity of 3 ab −1 are shown in Fig. 1 as solid and dashed lines, respectively, for m S = 200 GeV and three different DM masses, m χ = (5, 100, 300) GeV. 6 To obtain the projections, we used CheckMate with upscaled event numbers assuming that ATLAS measures the same distributions. Following [23] we further assume that the background error can be lowered by a factor of 4. Due to the nature of the process, radiating two DM particles from an internal mediator, interestingly the limits do not die off quickly 6 While with this choice the flavor model considered is fine, note that for m S 225 GeV strong bounds on the S-Yukawa couplings arise from the recent ATLAS search for resonant di-lepton production [22], which would exceed the projected limits of Fig. 1. Clearly, this can be avoided by moving either to the leptophobic or the hadrophobic scenario. when m χ > m S /2, allowing to test also this hierarchy of masses. As mentioned, further improvement could be reached by adjusting the analysis to the specific signature, e.g. by demanding two correlated jets in the final state. We leave the detailed study for future work.
We finally note that, although the final state looks similar to the one of Higgs to invisible searches in vectorboson fusion production, we found that the distribution of our signal in the main kinematic variables is very similar to the main backgrounds in that analysis and therefore no effective separation is possible there.
III. e + e − + E miss T

AT CLIC
An interesting proposal for a next high-energy e + e − collider facility to be built is the Compact Linear Collider (CLIC) at CERN. It would be the first mature realization of a collider with these characteristics and could start running in 2035. In the following, we will analyze the prospects to probe y S χ at the three foreseen stages of CLIC, stage I with √ s = 380 GeV, stage II with √ s = 1.5 TeV and stage III with √ s = 3 TeV. The corresponding luminosity goals are 1.0 ab −1 , 2.5 ab −1 , and 5 ab −1 , respectively [24,25].
To test the Z 2 -symmetric portal we propose a search in the e + e − + E miss T final state at CLIC, with the signal processes depicted in Fig. 2. The main irreducible background is [26] e + e − → e + e −ν ν , with the most important contribution coming from a ZZ intermediate state, while further backgrounds turn out to be negligible. For generating the signal and background samples at leading order, we employ again MadGraph5 aMC@NLO for the event generation, Pythia 8.1 for the hadronization and Delphes 3 for a fast detector simulation. The final analysis is performed with MadAnalysis 5 [27,28].
As it turns out, in the full flavor model, where S couples to electrons and quarks, the signal is very small for realistic couplings since the branching to quarks will strongly dominate (while simultaneously increasing significantly the total width). So we first focus on the hadrophobic case, with y S d = y S u = 0. 7 Still, we have to face a rather small signal with a sizable background, leading to weak constraints from a pure cutand-count analysis, in particular when the uncertainty in the background cross-section normalization is taken into account. Therefore we perform a shape analysis with a binned likelihood approach, making use of the fact that our signal has a peak-like structure in the m ll -variabledue to an on-shell S decaying to electrons -compared to a smoothly falling background. 8 To achieve a preliminary separation between signal and background, we apply the cuts given in Tab. I, where the 7 It would also be interesting to consider the di-jet final state at CLIC or to constrain the bi-quadratic portal at other colliders, however these analyses face their own challenges and will be left for future work. 8 In fact, the resonant diagram in the right panel of Fig. 2 largely dominates the cross section.
MET m e + e − pT (e) ∆R(e + e − ) θ(e + ) θ(e − ) m e + e − cut is applied to lower the impact of Z decays. In Fig. 3 examples of the shapes of signal and background after cuts and before fitting are shown for stage III. Here, the couplings y S e /Λ = 1.5/TeV and y S χ /Λ = 0.25/TeV are chosen to be close to the exclusion limit (see below).

A. Fitting Signal and Background
In order to use the m ll spectrum to discriminate signal and background, we first need sizable Monte-Carlo samples of both processes, where we generate 50.000 and 10 6 events, respectively. Since the signal shape depends on the width of S, it is simulated for various values of the latter, depending non-trivially on the input parameters (basically m S and y S e /Λ) given at the end of Section I. The resulting histograms are fitted to a 4-th order polynomial for the background and a simple Breit-Wigner distribution for the signal. Finally, the signal is characterized by the total number of events and the width of the Breit-Wigner distribution, allowing to easily test several couplings.

B. The Likelihood Function
To derive exclusion regions, we start with a binned Likelihood function [29] for the number of events n i , similar to the one used in CheckMate [20], and Here, S and B are the predicted numbers of signal and background events, respectively, while θ S,B are nuisance parameters incorporating the corresponding uncertainties ∆S and ∆B. Finally, the variation of the signal strength with the input parameters, given in Sec. I C, is parameterized by the signal-strength modifier µ, which is normalized for fixed y S e /Λ and fixed masses such that µ = (y S χ /Λ) 2 .
To test the compatibility of different values for the latter with data, we use the profile likelihood ratio [29] wereθ S (µ),θ B (µ) maximize L for the given value of µ, whileμ,θ S ,θ B correspond to the unconditional (global) maximum appearing in the denominator and are called unconditional Maximum Likelihood (ML) estimators. Here, the lower case accounts for the fact that we can only have a positive signal contribution.
Finally, for the numerical analysis it is convenient to use the test statistics [29] to set upper limits (with higher values corresponding to less compatibility), for which we use the python package iminuit [30].

C. P-Values
In the following we assume that the true underlying theory features µ = 0, i.e. we expect to see background only, and want to derive corresponding projected experimental exclusion regions on µ.
In general, to quantify the agreement between a (potentially) observed measurement and a signal hypothesis µ > 0, leading to a certainq µ,obs , the p−value is calculated, where f (q µ |µ ) is the probability density function (pdf) ofq µ under the assumption that the data is distributed according to a true µ = µ , while the subscript in the first argument denotes the hypothesis being tested. 9 As we want to derive the expected upper limits from future experiments, assuming no signal to be present, we will use the median value of the corresponding distribution, f (q µ |0), forq µ,obs . Finally, working at the 95% confidence level, we will solve for the value of µ that leads to p µ = 0.05.
To obtain the distributions f (q µ |µ ) without a large number of Monte Carlo simulations, we use the asymptotic formulas given in Ref. [29]. Those are valid for a sufficiently high number of events in each bin, which is fulfilled in our case. 10 While in the case µ = µ, f (q µ |µ) is given by a simple half-chi-square distribution, for obtaining the median ofq µ according to f (q µ |0) the socalled Asimov data set is used [29], where all estimators obtain their true values. This data set can be approximated via large MC simulations. Here we assume that our initial sets are large enough and use the fitted distributions as Asimov data. With this, the corresponding Likelihood-function and test statistics can be evaluated, which are denoted by L A and q µ,A . The variance, from which f (q µ |0) can be obtained, is then simply given by σ 2 A = µ 2 q µ,A , assuming background-only [29]. In practice we can however just use the Asimov value q µ,A for the median of [q µ |0], according to [29], and therefore the expected p−value for a signal hypothesis becomes with Φ the cumulative Gaussian distribution. In the end, p µ is evaluated for varying µ to find p µ = 0.05.

D. Resulting Limits
To establish the constraints on the model parameters, we have to translate the limits on µ into limits for the former. As mentioned before, for fixed y e S and thereby fixed width and shape of the m ee distribution, we have µ = (y S χ /Λ) 2 . For all limits we take a 5% uncertainty on the background normalization into account, i.e., σ B = 0.05 (while σ S is negligible).
In Fig. 4 we compare the reach of the three CLIC stages on the couplings, assuming m S = 200 GeV and m χ = 5 GeV. We observe that already at the first stage we would be sensitive to O(1/TeV) couplings, while at  the later stages the reach extends well beyond a TeV. In Fig. 5 the expected limits obtained for the same m S , but varying dark matter masses, are shown for CLIC stage II, which demonstrates that the sensitivity does not vanish for m χ /2 > m S .
We further note that direct searches for the mediator, e.g. in the e + e − final state, could break the degeneracy between the two couplings. It might well happen that the mediator would first be found via such a search, however then the present analysis would be crucial to investigate the structure of the dark sector.

IV. DARK MATTER PHENOMENOLOGY
For m χ m S , the DM relic density is set via the processχχ → SS, while for smaller dark matter masses it is always far above the measured value since no decay channel is kinematically allowed (the s-channel decay induced by v S > 0 is found to be negligible, even in the resonance region). The viable parameter region, featuring 0.11 < h 2 Ω DM < 0.13, is shown as a blue band in Fig. 6 in the m S −m χ plane, where we set y S χ = 2.25. Light mediators m S < 200 GeV, below the green line, are already excluded by XENON1t [31] and heavier once will be tested in future experiments like LZ [32] (red line) and DAR-WIN [33] (remaining region). The dominant contribution to direct detection rates arises from tree-level s-channel exchange of S with the up and down quarks and therefore vanishes in the hadrophobic case. Since v S ∝ 1/y S f , the cross section is independent of the Yukawa couplings. All numerical results have been obtained with micrOmegas 5.0.8 [34].
Finally, the required y S χ in dependence on m χ is shown in Fig. 7 for m S = 200 GeV. Note that also the relic density is independent of the values of y S u (or y S e ), which do not enter the dominant annihilation amplitude. We find that, unless the electron S−Yukawa coupling is very small, most of the viable parameter space will be tested at CLIC.