Signatures from Scalar Dark Matter with a Vector-like Quark Mediator

We present a comprehensive study of a model where the dark matter is composed of a singlet real scalar that couples to the Standard Model predominantly via a Yukawa interaction with a light quark and a colored vector-like fermion. A distinctive feature of this scenario is that thermal freeze-out in the early universe may be driven by annihilation both into gluon pairs at one-loop ($gg$) and by virtual internal Bremsstrahlung of a gluon ($q \bar{q} g$). Such a dark matter candidate may also be tested through direct and indirect detection and at the LHC; viable candidates have either a mass nearly degenerate with that of the fermionic mediator or a mass above about 2 TeV.


I. INTRODUCTION
In recent years, dark matter in the form of Weakly Interacting Massive Particles (WIMPs) has become the leading particle physics candidate. The most salient feature of WIMP dark matter is the prediction of a relic abundance via thermal freeze-out that is in rough agreement with the measured value, when the WIMP annihilation cross section into Standard Model (SM) particles is of the order of 1 pb, thus suggesting parameters in the dark sector (WIMP mass and couplings) comparable to those in the electroweak sector. Furthermore, many particle physics models addressing the electroweak hierarchy problem contain, among many new states, a dark matter candidate with the characteristics of a WIMP. While the lack of experimental evidence for these extra states is disfavoring some of these schemes, the freeze-out mechanism still stands as one of the most natural and elegant mechanisms to explain the origin of the dark matter.
An exciting feature of the WIMP dark matter paradigm, and arguably one of the reasons for its popularity, is the possibility of detecting experimental signals of processes induced by the weak coupling of the dark matter with the Standard Model particles. Numerous experiments are currently searching for dark matter particles, either through direct and indirect detection or at colliders, and are already providing fairly stringent limits on the rates of the processes relevant for each search strategy. Unfortunately, the limits on the fundamental parameters of the model which can be derived from the null searches, as well as the complementarity among the various search strategies, are highly model dependent. On the other hand, the main features of the dark matter phenomenology of a given model can be, in many instances, captured by considering only a subset of the new fields and new parameters, namely by considering "simplified models".
The simplest among all the simplified WIMP models is to extend the Standard Model with a singlet real scalar, S, and with a discrete Z 2 symmetry, unbroken in the electroweak vacuum, under which the singlet scalar is odd while all the Standard Model fields are even [1][2][3][4], and which ensures the dark matter stability. 1 The singlet scalar interacts with the Standard Model 1 Other parity assignments are possible. For instance stability of scalar DM may be explained through embedding in SO(10) using matter parity P = (−) 3(B−L) , which is a remnant of a gauge symmetry [5]. In this framework, SM fermions are as usual in a 16 SM (odd under P ) while the SM Higgs is in 10 H (even). If P is unbroken, a possibility is to embed scalar DM in a 16 DM , which contains a singlet scalar and a scalar doublet (aka inert doublet). Either one may play the role of DM [5]. Notice that Yukawas of the form of Eq. (2) may be obtained by coupling the 16 DM to 16 SM through fermions in the 10 → 5 ⊕5, which contains a SU (2) L singlet vector-like quark, akin to the d R , as well as a leptonic SU (2) L doublet. Of course more work would particles only through a renormalizable quartic coupling to the Brout-Englert-Higgs doublet (Higgs for short) thus the model only contains one new field and two new parameters, m S and λ; the phenomenology of this very simple scenario has been studied in many works (see for instance [12][13][14][15][16][17][18]).
In this paper we will investigate the impact on the dark matter phenomenology of extending this model by one vector-like fermionic field, ψ, also odd under the discrete Z 2 symmetry, and which couples to the singlet scalar S and a right-handed SM fermion f R via a Yukawa interaction (the discussion for a left-handed SM fermion is completely analogous). The Lagrangian of the Z 2 -odd fermionic sector reads: The gauge invariance of the Yukawa term implies that the fermion must be a SU (2) L singlet and have the same hypercharge as the SM fermion it couples to. Besides, the vector-like nature of the fermion ensures the cancellation of the gauge anomalies. 2 From a more fundamental perspective, notice that this model is also a simplified version of a scenario with large extra dimensions (see [19] and references therein).
As we will argue, when the strength of the "vector-like portal" interaction is larger than the strength of the "Higgs-portal" interaction, the phenomenology of the model presents quite distinctive characteristics which make this model rather unique. The phenomenology of the vector-like portal interaction has been partially discussed in, e.g., [19][20][21][22][23][24][25][26][27]. The case in which f R is a light lepton was first discussed in [23,24]. There, it was shown that the annihilation of S through the portal of Eq. (2) into two SM fermions is dominated, in the limit m f m S , by the d-wave, thus σv ∝ v 4 , and that the virtual internal Bremsstrahlung (VIB) in the annihilation may lead to rather intense gamma ray spectral features (see also [28][29][30][31][32] for relevant literature on VIB in an analogous framework with Majorana dark matter). Further phenomenological be necessary to make this a phenomenologically viable framework. For recent works along this direction, see [6][7][8][9][10][11]. 2 Taking ψ chiral instead requires introducing more fields, but is otherwise an obvious extension (with a possible caveat regarding the mass of the ψ, which may have to come from couplings to the Higgs and thus may not be arbitrarily large). aspects of this "leptophilic" scenario were studied in [26] and [27]. In particular [26] analyzed constraints both from indirect detection through gamma rays and from colliders.
In the present work, we will study the possibility that the dark matter couples to a light quark, assuming for simplicity minimal flavour violation. An interesting feature of this scenario is that, due to the d-wave suppression of the annihilation rate into a quark-antiquark pair, the relic abundance is driven either by the annihilation of S into two gluons 3 (gg) and into VIB of a gluon (qqg), or by the co-annihilation of the colored ψ particle. The latter process can be significantly affected by the exchange of gluons between the non-relativistic fermion mediators, a phenomenon that may lead to Sommerfeld enhancement of the annihilation crosssection. The correct description of the phenomenology of the model then requires the inclusion of higher-order effects, as well as non-perturbative effects.
The fact that the fermionic mediator ψ is colored also opens interesting possibilities for direct detection and collider searches 4 . Specifically, we consider the constraints set by the LUX [35] experiment and the prospects for detection by XENON1T [36]. To this end, we take into account the effective coupling at one-loop of the S to gluons, recently calculated in [37]. Furthermore, the model can be tested at colliders through the production of the mediator particles ψ, which are both colored and electrically charged, and which subsequently decay producing a signature of two or more jets plus missing transverse energy E T . Altogether, direct detection and production of the fermionic mediator at colliders put the strongest constraints on our scenario. Finally we also consider indirect searches, in particular constraints on continuum gamma-ray signals from dwarf spheroidal galaxies (dSphs) and gamma-ray spectral features, based on data from Fermi-LAT and H.E.S.S., and on anti-protons in cosmic rays, using the PAMELA data. While overall less constraining, these data sets will allow us to close a narrow region of the parameter space that otherwise would be left open by current direct and collider searches. Similar analyses applied to other dark matter scenarios have been pursued in 3 The possibility that the annihilation into gluon pairs could drive the dark matter freeze-out was previously considered in [33], using a model-independent effective operator framework. The present model provides an explicit UV realization of this scenario. 4 As the couplings we consider violate flavour symmetry, it is legitimate to ask whether our scenario may be made consistent with constraints on flavour changing neutral processes. This issue has been addressed in [34] where two mechanisms to suppress such processes (degeneracy and alignment) are discussed. As these mechanisms directly apply to our model, we do not repeat the argument but refer to [34] for more details.
This article is organized as follows. We first determine in Sec II the dark matter relic abundance, taking into account higher-order processes involving gluons, as well as the Sommerfeld enhancement for co-annihilation of the fermionic mediator. Then, we study in Secs. III, IV and V the constraints from direct searches, collider searches and indirect searches, respectively.
We finally draw our conclusions in Sec. VI.

II. DARK MATTER RELIC DENSITY
As will be shown below, the relic dark matter abundance is determined either by the one-loop annihilation into gluon pairs (gg) or the VIB in the annihilation (qqg), or by the co-annihilation with the colored partner. We first discuss the perturbative aspects of the annihilation, and then the Sommerfeld effects in the co-annihilation processes associated to the exchange of gluons between particles in initial states.

A. Perturbative analysis
We will focus on the scenario where the vector-like portal interaction dominates over the Higgs portal interaction, namely we assume that the quartic coupling is small enough to play no role. In this case, the lowest order dark matter annihilation channel is SS → qq, with the annihilation cross-section given by σv where we have kept only the leading terms in an expansion in m q m S and relative velocity v. Here, r ≡ m ψ /m S > 1 denotes the mass ratio between the vector-like mediator and the dark matter particle. It follows that in the limit m q → 0, which is relevant for the case of dark matter coupling to a light quark, the annihilation of dark matter into a pair of quarks is d-wave suppressed [23,24]. While the absence of the s-wave component in the limit m q → 0 is analogous to the well-known helicity suppression in the case of annihilation of a Majorana fermion pair into light fermions, the vanishing of the p-wave component is unique to the annihilation of two real scalar particles; see [24] for a more detailed discussion. Due to the strong suppression of the tree-level two-to-two cross-section, it is crucial to take into account higher-order effects in the calculation of the relic density. There are two relevant classes of higher-order processes which contribute to the total annihilation cross section: SS → qqV , i.e. the internal Bremsstrahlung of a gauge boson V , with V = γ, Z, g, as well as SS → V V , the one-loop annihilation of dark matter into a pair of gauge bosons V, V . Representative diagrams for these higher-order annihilation channels are shown in Fig. 1, together with the treelevel annihilation process SS → qq. In contrast to Eq. (3), the higher-order processes feature a non-vanishing s-wave component in the limit m q → 0, at the price of being suppressed by an additional gauge coupling and reduced phase space in the case of internal Bremsstrahlung, or by a loop suppression factor for annihilation into gauge bosons. Analytical expressions of the corresponding cross sections can be found in [23,24,26,27].
To calculate the relic density, we employ the micrOMEGAS.4 package [44], which we modified in order to include the additional annihilation channels. 5 We find that in the whole parameter space either the internal Bremsstrahlung processes or the one-loop annihilations give a larger contribution to the total annihilation cross section at freeze-out than the annihilation into qq given by Eq. (3). Due to the different scaling with the mass ratio r of the two higherorder processes, namely σv ∝ 1/r 8 for internal Bremsstrahlung and σv ∝ 1/r 4 for one-loop annihilations (see [23,24,26,27]), the former dominates the total annihilation cross section for small mass ratios (concretely, r 2.5), while the latter is more important for large mass ratios.
Furthermore, if the fermionic mediator is close in mass to the dark matter particle, coannihilation processes, i.e. annihilations of S and ψ, ψ and ψ, or ψ andψ [45], can play a significant role in the calculation of the relic density. As is well known, the co-annihilation channels have rates exponentially suppressed when ∆m ≡ m ψ − m S T f.o. , yet, for small mass splittings, they can dominate over the self-annihilation channels, which have in this model suppressed rates. We employ micrOMEGAS.4 in order to include all relevant co-annihilation processes in our numerical calculation, and we find that co-annihilations give a contribution to Ω DM of at least 5% for mass splittings r 1.5, 1.35, 1.25 at m S = 10 GeV, 100 GeV, 1 TeV, respectively. It is important to note that some of the co-annihilation processes as e.g. ψψ → gg are purely given by gauge interactions; in particular, they are independent of the Yukawa coupling y. Consequently, if m S and r are small enough, these annihilation channels can suppress the relic density below the observed value, irrespectively of the value of y. 6 For the parts of the parameter space where co-annihilations are not relevant, the relic density is set by the higher-order processes introduced above, and consequently rather large Yukawa couplings are typically necessary for matching the observed relic density. As it will be discussed in the rest of this work, this implies potentially observable signals in direct and indirect detection, as well as in the production of the mediator at colliders.

B. Sommerfeld corrections to the co-annihilation processes
The annihilation of two colored particles (for our purposes, with initial states ψψ,ψψ or ψψ) can be significantly affected by the non-perturbative Sommerfeld effect, induced by the multiple exchange of gluons in the initial state. For the quantitative treatment of the Sommerfeld effect, we closely follow the method presented in [43,46], which we briefly recapitulate here.
Expanding the perturbative annihilation cross section as the Sommerfeld corrected cross section can be written as Here, S l is the Sommerfeld factor associated with the partial wave l, which takes the form [47,48] In these expressions, v is the relative velocity between the two annihilating particles, and the coupling α parametrizes the strength of the effective QCD potential between the annihilating particles. For two annihilating particles, with representations R and R under SU (3) c , the non-abelian matrix potential between the particles can be diagonalized by decomposing the direct product R ⊗ R = Q Q as a sum of irreducible representations Q. This gives rise to the effective potential (see e.g. [46]) where we evaluate α s at an energy scale equal to the momenta of the annihilating particles, p = mv/2 and the C i are the quadratic Casimir operators associated with the representation i.
For the model under consideration, we will be interested in processes involving the annihilation of colored fermion triplets giving rise to 3 ⊗3 = 1 ⊕ 8 and 3 ⊗ 3 = 6 ⊕3. The associated Casimir operators are given by C 1 = 0, C 3 = C3 = 4/3, C 6 = 10/3 and C 8 = 3. The effective QCD potentials then read and are hence attractive for the singlet and the anti-triplet two particles states while they are repulsive for the others. In the following, we will refer to S l , S As the Sommerfeld enhancement depends on the initial color state, one has to determine the relative probabilities of annihilation in a given color state, separately for every annihilation channel. As discussed in [46], using tensor decomposition, one can show that ψψ → gg occurs with probability 2/7 (5/7) through a singlet (octet) state, giving rise to the total Sommerfeld For the annihilation channel ψψ → qq due to the Yukawa interaction of Eq. (2), a similar calculation yields On the other hand, the annihilations of ψψ into γg, Zg (γγ, γZ, ZZ, W W, Zh) correspond to a pure octet (singlet) initial state, due to color conservation. The associated Sommerfeld factors are thus given by Finally, we also consider the co-annihilation process ψψ →qq. When the dark matter has a tree-level Yukawa interaction with the quark q, the co-annihilation process is mediated by the t-channel exchange of a fermionic mediator, as well as by the s-channel exchange of a gluon.
If this is not the case, only the s-channel process is relevant. As a consequence, the SU (3) c representation of the initial state crucially depends on the concrete quark produced in the final state. Namely, if only the s-channel mediates the co-annihilation, the initial state is a pure octet. On the other hand, if the s-and t-channels interfere, the initial state is combination of singlet and octet representations, therefore the total Sommerfeld factor can not be expressed in a simple form analogously to Eqs. (9)-(11); instead, we fully decompose the squared matrix element into the part corresponding to a singlet and an octet initial state, and multiply each contribution with S (1) l and S (8) l , respectively.

C. Results of the relic density calculation
In this work, we consider the case of dark matter coupling either to the right-handed up-quark u R , or to the right-handed down-quark d R . Both for the annihilation and the co-annihilation processes, the dominant contribution to the annihilation cross section arises either from the Yukawa interaction or from the QCD interaction, the latter being identical for coupling to u R and d R . Hence, up to sub-percent corrections arising from diagrams involving electromagnetic or weak couplings, the Yukawa coupling corresponding to the observed relic density is identical in both of these scenarios. Note however that we fully include all sub-leading effects in our numerical calculations.
In Fig. 2 we show the viable parameter space in the plane spanned by the dark matter mass m S and the relative mass difference r − 1, as well as in the plane spanned by m S and the Yukawa coupling y. The color gradient corresponds to the value of the Yukawa coupling in the left panel and to the value of r − 1 in the right panel required to match the observed relic abundance Ω DM h 2 0.12. In both panels, the region at small m S and/or small r, i.e. the dark grey in the left panel, corresponds to parameters for which the co-annihilation processes involving only gauge interactions are sufficient to suppress the relic density below the observed value. Besides, in the left panel, the light grey region at large m S and/or r corresponds to Yukawa couplings larger than y = 6, which we choose as perturbativity limit. 8 Notice that for the most degenerate scenario discussed in this work, r = 1.01, the Sommerfeld corrections lead to a change in the relic density of up to 15%. It is also worth mentioning that the Sommerfeld effect can both increase as well as decrease the relic density, depending on the dark matter mass and the mass splitting. This is due to the different sign of the effective coupling α depending on the initial color state of the annihilation process (c.f. Section II B).

III. DIRECT DETECTION CONSTRAINTS
The interaction of dark matter particles with nucleons leads to potentially observable signatures in direct detection experiments. In the model under study in this paper, this interaction is described by the Feynman diagrams shown in Fig. 3. The diagrams in the upper row correspond to the tree-level scattering off a light quark q, being in our case q = u or d, through the exchange of a fermionic mediator, while the one-loop diagrams depicted in the lower row lead to an effective interaction of the dark matter particle with gluons.
The effective interaction of S with a light quark q can be cast as the sum of a scalar and a twist-2 contribution [49]: By integrating out the mediator ψ in the diagrams shown in the upper row of Fig. 3, and matching the result to the effective Lagrangian in Eq. (12), we obtain where q = u (q = d) for dark matter coupling to u R (d R ).
On the other hand, the diagrams in the lower row of Fig. 3 induce an effective coupling of dark matter to gluons, which has been recently discussed in [37]. Similarly to the well-studied simplified model of Majorana dark matter with a colored scalar mediator [50], the scattering amplitude can be decomposed into a short-distance and a long-distance contribution. This separation is made according to the momentum scale dominating the loop integration: the part of the amplitude arising from loop-momenta of the order of the mass scale of the heavy particles (i.e. the dark matter particle or the fermionic mediator) leads to the short-distance contribution, while the part arising from loop-momenta of the order of the light quark masses, leads to the long-distance contribution. The latter involves values of the strong coupling constant at a non-perturbatively small momentum scale, and hence can not be reliably calculated in perturbation theory. However, this contribution is implicitly contained in the parton distribution function of the light quarks in the nucleon, and hence only the short-distance contribution must be taken into account in the computation of the one-loop diagrams. Defining the effective Lagrangian for the dark matter-gluon interaction as the short-distance contribution to C g S reads, assuming the limit m q m ψ − m S [37], From the effective Lagrangians to the partons, Eqs. (12) and (15), one can calculate the effective spin-independent coupling of the dark matter particle S to a nucleon N , the result being [49] f where f (N ) T G are mass fractions and q (N ) (2),q (N ) (2) are the second moments of the parton distribution functions; for our numerical analysis we use the values in [37]. We show in Fig. 4, left panel, the effective dark matter-proton coupling f p /m p , expressed in units of y 4 /m 2 S (which is the common pre-factor to each of the terms in Eq. (17)), as a function of r − 1, for the case of dark matter coupling to u R . As apparent from the plot, the coupling of dark matter to gluons, corresponding to the last term in Eq. (17), interferes destructively with the contributions arising from the dark matter scattering off quarks, leading to a vanishing dark matter-proton coupling at r 3.0. A similar behavior is found for f n , although the total destructive interference occurs at a slightly different mass splitting r 2.3.
Finally, from the effective dark matter-nucleon couplings f p and f n , one can calculate the total spin-independent cross section for dark matter scattering off a nucleus with charge Z and mass number A. At zero momentum transfer, which is the regime relevant for direct detection experiments, it is given by the coherent sum of f p and f n : We find that, among all current experiments, the LUX experiment [35] gives, in most of the parameter space, the best limits to the simplified model under scrutiny in this paper 9 . LUX 9 After this work has been completed, the LUX collaboration presented a reanalysis of the 2013 data, leading to improved upper limits on the dark matter scattering cross section [51]. However, none of our conclusions would change qualitatively when including the updated limits.
is based on a xenon target and, as is well known, has least sensitivity when f n /f p −0.7, the precise value being dependent on the xenon isotope considered. In this case, commonly referred to as "maximal isospin violation", the total scattering cross section vanishes. We show in the right panel of Fig. 4 the ratio f n /f p as a function of the mass splitting, for the case of dark matter coupling to u R . As follows from the plot, in this scenario maximal isospin violation occurs at r 2.6 for dark matter coupling to u R (this occurs at r 3.3 for dark matter coupling to d R ). Hence, we fully take into account the actual value of f n /f p predicted by the model by defining the effective dark matter-proton cross section σ eff p via where ξ i is the natural relative abundance of the xenon isotope i. In that way, σ eff p can be interpreted as the dark matter-proton scattering cross section that under the assumption of f p = f n would lead to the same number of events in LUX as the actual dark matter-proton cross section σ p , when taking into account the ratio f n /f p predicted by the model. Consequently, σ eff p is the quantity that can be compared to the upper limit derived by the LUX collaboration [35], which was derived under the assumption that f p = f n .
In Fig. 5, we show the effective dark matter-proton scattering cross section σ eff p in two different projections of the thermal parameter space. In the left panel, the color gradient corresponds to different values of the relative mass difference between the vector-like mediator and the dark matter particle. Let us emphasize that for every dark matter mass and mass splitting, we fix the Yukawa coupling y by the requirement of matching the relic abundance, cf. section II C. As it can be seen from the plot, there is a huge scatter in the possible values of the scattering cross section over many orders of magnitude. Generally speaking, for a fixed dark matter mass, the scattering cross section gets larger for more and more degenerate scenarios, i.e. smaller r, due to the resonant enhancement of the dark matter-nucleon coupling f N in the limit r → 1, which follows from Eqs. (14), (16) and (17). However, as it can be inferred from to a fixed value in the limit of r 1. This can easily be understood by noting that in this limit, the most relevant process for the freeze-out of dark matter is the one-loop annihilation of dark matter into gluons, scaling as 1/r 4 , which is precisely the same asymptotic behaviour of the scattering cross section for r 1, as it can be seen from Eqs. (14), (16), (17) and (18).
In the left panel of Fig. 5, we also show the current upper limit from LUX, as well as the projected upper limit for the XENON1T experiment, taken from [36]. In addition, we also show in the plot the ultimate reach of (non-directional) direct detection experiments, given by the scattering cross section corresponding to coherent neutrino scattering in the detector [52]. As apparent from the plot, for m S 200 − 300 GeV, large parts of the parameter space are already excluded by LUX, and for the most degenerate scenarios even dark matter masses around 2 TeV are already ruled out. In the near future, XENON1T will continue closing in the region of the parameter space of thermal dark matter, in particular for dark masses above the TeV scale, which will be difficult to probe at the LHC (see the discussion in the next chapter).
Finally, in the right panel of Fig. 5 we show the effective dark matter-proton scattering cross section as a function of r − 1, with the color coding corresponding to the Yukawa coupling y, and fixing the dark matter mass by the relic density requirement. In this projection of the parameter space, the suppression of the effective cross section at r 2.6 becomes apparent, due to the maximal isospin violation occurring at this value of the mass ratio.

IV. CONSTRAINTS FROM SEARCHES AT THE LHC
The model discussed in this work can be efficiently probed at the LHC through the production of the fermionic mediator ψ. After being produced, the mediator decays promptly into the dark matter particle S and a light quark; hence, the signature of interest consists of (at least) two jets plus missing transverse E T . In the following, we first discuss the relevant contributions to the production cross section of the fermionic mediator ψ, and subsequently we derive constraints on the model by reinterpreting a recent ATLAS search for multiple jets plus missing E T .

A. Production of mediator pairs
We consider the production of mediator pairs ψψ, ψψ, andψψ in proton-proton collisions at a center-of-mass energy √ s = 8 TeV. 10 . The relevant Feynman diagrams, shown in Fig. 6, can be divided into two categories: first, a ψψ pair can be produced from a qq or gg initial state by QCD interactions, as shown in the upper row of Fig. 6. Secondly, for the case of dark matter coupling to u R , a pair of ψψ, ψψ orψψ can be produced from a uu, uū, orūū initial state, respectively, through the t-channel exchange of the dark matter particle S, as shown in the lower row of Fig. 6. An analogous statement holds for dark matter coupling to d R , by exchanging u with d in the initial states. The amplitudes corresponding to the QCD induced diagrams are independent of the dark matter mass and the Yukawa coupling y, in contrast to 10 In addition to the processes pp → XX, with X = ψ,ψ considered here, the production of mediators through the process pp → XS provide extra contributions to the multiple jets plus missing E T analysis. Assuming y = y thermal , we have checked that the total production cross sections satisfy σ(pp → XS) < σ(pp → XX) and that the enhancement in the production cross-section obtained when including σ(pp → XS) does not affect significantly the limits that are presented in Fig. 8. the processes arising from the Yukawa interaction.
In the model discussed in this work, the requirement of correctly reproducing the observed dark matter abundance via thermal freeze-out can involve rather large Yukawa couplings, as can be seen in e.g. Fig. 2. On the other hand, the strong coupling g s 1, therefore the production of mediator pairs can be dominated either by the processes involving the Yukawa coupling or by the strong coupling. Furthermore, for dark matter coupling to u R (d R ), the production channel uu → ψψ (dd → ψψ), which is driven solely by the Yukawa coupling, can have a larger cross section compared to processes involving a qq or gg pair in the initial state, due to the large parton distribution function of the up-quark (down-quark) in the proton. In our analysis, we calculate the leading-order cross section for the final states ψψ, ψψ,ψψ using CalcHEP [53], selecting the cteq6l parton distribution function. In order to include next-to-leading order corrections, we parametrize the total cross section as The K-factor has not been fully calculated for this model, consequently we treat K as a source of uncertainty, varying it in the range [0.5, 2]. 11 11 For comparison, the K-factor for tt production, which is analogous to the QCD production of ψψ, is 1.7 [54]. We show in Fig. 7 the production cross section for the various initial states for the exemplary choice of dark matter coupling to u R , and choosing m ψ = 500 GeV as well as K = 1. The red curves were calculated for the value of the Yukawa coupling that leads to the observed dark matter abundance via thermal freeze-out, y = y thermal , while the blue dashed curve shows, for comparison, the predicted cross section when taking into account only pure QCD processes in the production of mediator pairs. Notice the rise of the total cross section for large values of the mass ratio, that is entirely due to the Yukawa induced processes. In particular, for m ψ /m S 1.2, the production channel uu → ψψ dominates, due to the large thermal Yukawa coupling as well as due to the large parton distribution function of the up-quark in the proton, as discussed above. From Fig. 7, it also follows that the cross section for the process uū → ψψ is lower than the pure QCD cross section in the range 1.1 m ψ /m S 1.3. This is due to the destructive interference of the QCD and Yukawa driven contributions to the process.

B. Constraints derived from the multijet ATLAS analysis
The production of ψ can be detected at the LHC through its decay into a quark and the dark matter particle S. Ignoring the higher-order corrections, the event topology always consists of two jets arising from the decay of the two mediators, together with missing E T associated to the dark matter particles (this process was considered in [55], neglecting, however, the impact of co-annihilations and higher-order effects in the calculation of the dark matter abundance). On the other hand, the radiation of quarks or gluons from the initial, final or an intermediate state can lead to additional jets which could also be detected. Considering these extra jets from higher-order processes is important due to two reasons: first, if the absolute mass splitting m ψ − m S is below ∼ 50 − 100 GeV, the jet arising from the decay of ψ is too soft to be detected at ATLAS or CMS. In that case, the additional emission of one or more hard jets is necessary for the detection of the event (for an analysis of the constraints of the model employing monojet searches, see [41,55]). Secondly, even if the mass splitting of the dark matter particle and the mediator is sufficiently large, the signal-to-background ratio of the n-jet + / E T signal can be larger for n > 2 compared to n = 2, depending on the characteristics of the relevant Standard Model background processes. Consequently, even if the absolute cross section is smaller, the constraints derived from the multijet processes can be more stringent than those obtained from the lowest order two-jet topology.
In order to confront our model to available LHC data, we make use of the ATLAS search in each signal region, which can be found in [56], are denoted as S 95 obs (S 95 exp ). On the other hand, for a given signal region i, the number of expected signal events can be cast as where σ is the total production cross section for pairs of ψ as discussed in section IV A, and L = 20.3 fb −1 is the luminosity. Moreover, the efficiency i is the probability that a given event passes all the selection cuts corresponding to the signal region i. This quantity depends on the topology of the involved Feynman diagrams, and hence has to be recalculated for a given model. For that purpose, we have implemented our model in FeynRules [57], and use MadGraph [58] in order to simulate events for the production of mediator pairs with up to two additional jets. The partonic events are fed to PYTHIA 6 [59], which simulates the showering and hadronization process. As usual in this context, in that step one has to take care of the possible double counting in the additional jets, which should either be included at the matrix element level in the partonic event, or as part of the showering process. In order to ensure a smooth transition between the two regimes, we employ the MLM matching scheme, using the parameters Qcut = m ψ /4 and SHOWERKT=T. We explicitly checked that with this choice the differential jet distributions are smooth, as required for a physical meaningful matching scheme.
Finally, in order to apply the cuts corresponding to the ATLAS search, we use CheckMATE [60], which also implements a detector simulation. Concretely, for given values of m S , m ψ , and the corresponding Yukawa coupling y = y thermal (m S , m ψ ), we simulate N ev events as described above, and for every signal region i CheckMATE determines the number of events after cuts passing all the cuts. Then, the efficiency can be obtained through i = N (i) after cuts /N ev . As only a finite number of events can be generated, the efficiency obtained in this way is subject to Monte Carlo uncertainties. In order to take these into account in a conservative way, we replace S i by the 95% C.L. lower limit on S i , given by S i − 1.96∆S i , where ∆S i is the 1σ statistical error on the number of expected signal events 12 . Then, the 95% C.L. upper limit on the production cross section following from the observed (expected) number of events in the signal region i is given by Depending on the specific point in the parameter space, a different signal region provides the most stringent constraint. In order not to bias the final result by making a choice of the optimal signal region based on an over-or under-fluctuation, we choose for a given m S and m ψ the signal region giving rise to the smallest expected upper limit σ exp,i , and then use the observed upper limit σ (95) obs,i in that signal region as our final upper limit on the production cross section. As an example, this upper limit is shown in Fig. 7 for m ψ = 500 GeV, together with the production cross section predicted by the model. Our results for the case of dark matter coupling to u R show that the ATLAS data rules out dark matter masses below ∼ 800 GeV for r 1.5, upon imposing our perturbativity condition y < 6. This conclusion, furthermore, is rather insensitive to the concrete choice of the K factor. The ATLAS constraints for the scenario of coupling to d R are slightly weaker, due to the smaller parton distribution function of the down-quark in the proton, compared to the up-quark. We also show for completeness the limits from production of vector-like quarks at the Z resonance at LEPI (yellow dotted) and at LEPII (yellow solid).
For the former, we imposed m ψ m Z /2. For the latter, we assumed that the constraints are essentially kinematical and amount to m ψ 100 GeV, similar to LEPII results on squarks searches [61]. Besides, we also include in the plots, as a light red region enclosed by a solid for small values of r. A future 100 TeV proton-proton collider would perfectly complement the search for large r, and in combination with XENON1T, should be able to close in on the parameter space of the model for dark matter masses as heavy as a few tens of TeV.

V. INDIRECT DETECTION CONSTRAINTS
Indirect dark matter searches, using in particular gamma-rays and antiprotons as messengers, provide a complementary avenue to probe our model.
The model under consideration produces gamma-ray signals in the form of sharp spectral features or in the form of a continuum. The sharp spectral features arise mainly from the annihilation final states uūγ and γγ (the contribution to the total spectrum from γZ is a factor tan 2 θ W 0.30 smaller and will be neglected here). The photon multiplicities are analogous to those produced when the dark matter couples to a lepton, which were discussed in [23,24,26,27]. The final state uūγ gives the dominant contribution to the energy spectrum when the scalar DM and the fermion mediator are very close in mass, while γγ is more important in the opposite case; for intermediate regimes, the spectrum will be a linear combination of these two, weighted by the corresponding branching fractions. Therefore, for a given dark matter mass, the observational upper limits on uūγ and γγ bracket the upper limit on σ uūγ + 2σ γγ for any value of the fermion mediator mass. The total annihilation section into sharp gamma-ray spectral features, σ uūγ + 2σ γγ , is shown in Fig. 10 On the other hand, the continuum of gamma-rays is mainly produced by the final states uūg and gg (we neglect in our analysis the contributions from uūγ and γZ). Similarly to the discussion above, for small mass ratios the energy spectrum of continuum photons is dominated by annihilations into uūg, while for large mass ratios, by annihilations into gg; the corresponding spectra from annihilation in these two extreme cases are shown in Fig. 9. 13 Again, for intermediate regimes the photon spectrum is a combination of these two, and therefore one can bracket the upper limit on the cross section for any mediator mass by considering the upper limits on the channels uūg and gg. The total annihilation cross section for these two processes, σ uūg + σ gg , is shown in Fig. 10, right panel, together with the band bracketing the upper limits on σū ug + σ gg from the Fermi-LAT observations of dSphs [68], and which is bounded by the upper limits on σū ug and σ gg . 14 As can be seen from the Figure, the Fermi-LAT limits on continuum gamma-ray emission from dSphs put constraints on the parameter space of the model for m S 150 GeV. 13 We also show in the same figure the gamma-ray spectrum from the uū final state. The VIB differential cross section is peaked toward E g ∼ E u/ū ∼ m S at the parton level. This explains why the spectrum in gammarays and anti-protons is roughly half-way between that from gg and uū in Figs. 9 (see also Ref. [67]). This is illustrated for a compressed mass spectrum (r = 1.01) but we have checked that this feature also holds for larger values of r. 14 The Fermi-LAT collaboration has not published limits on these two final states, therefore, we estimate these limits following the methodology pursued in [32], using the 95 % CL exclusion limits on DM annihilation into qq from [68]. Concretely, we use the fact that the total gamma-ray flux, e.g. for the uū channel, is proportional to σv uū N uū γ , where N uū γ is the number of gamma rays per annihilation within the energy range of Fermi-LAT, 0.5 GeV< E γ < 500 GeV. We then rescale the Fermi-LAT limits on σv uū to get σv gg,uūg = σv uū N uū γ /N gg,uūg γ , where we have determined the number of gamma-rays N uū,gg,uūg γ using PYTHIA 8.1 [69]. The annihilation channels uūg and gg also produce an antiproton flux, which is currently constrained by the PAMELA data [70] (preliminary data from AMS-02 have been presented in [71]). We show in Fig. 10, right panel, the band bracketing the upper limit on σū ug + σ gg using the limits on σū ug derived in [31,72] and the limits on σḡ g derived in [33]. It follows from the should be borne in mind, though, that the antiproton limits are rather sensitive to uncertainties related to cosmic-ray propagation (these uncertainty may be reduced when more cosmic ray data will be released by the AMS-02 collaboration).
Finally, we confront in Fig. 8

VI. CONCLUSION
In this paper we have studied a simplified dark matter scenario in which a real scalar dark matter candidate S couples to light standard model quarks through a Yukawa interaction involving a new vector-like fermion mediator ψ. In order to compute the relic dark matter abundance via thermal freeze-out, we have taken into account annihilation into three body final states involving internal Bremsstrahlung as well as loop processes. We have also accounted for Sommerfeld corrections of the co-annihilation processes. The latter give rise up to a 15% change in the relic density (enhancing or decreasing the one calculated at perturbative level). In addition, a distinctive feature of our dark matter scenario is that annihilations into gg and qqg can be the dominant contributions to the dark matter annihilations both in the early universe and today.
We have investigated three complementary roads in order to test our dark matter scenarios: direct detection dark matter searches, collider searches and indirect detection through annihilation into antiprotons and gamma rays. Current constraints from the LUX experiment exclude most of the viable parameter space below masses of a hundreds of GeV (as well as a small island in the TeV range). The reach of the future XENON1T experiment extends to 10 TeV for low values of m ψ /m S − 1. In addition, we have analyzed how the production of the fermionic mediator at colliders could complement such constraints. To this end, we made use of an AT-LAS search for 2-6 jets + missing E T using pp collision data at √ s = 8 TeV with L = 20.3 fb. In order to compute the efficiencies corresponding to our particle physics model we have simulated events for the production of mediator pairs using Madgraph and applied the cuts corresponding to the ATLAS search using CheckMATE. Our results are the following: for masses between ∼ 30 (100) GeV and 1 TeV, such collider searches allow us to probe viable dark matter scenarios that are beyond the reach of current (future) direct detection searches. Combining the multijet + missing E T ATLAS search limit with LUX constraints, real scalar dark matter coupling to u or d quarks is excluded for masses in the range 5−300 GeV. Moreover, the ATLAS limit extends up to m S ∼ 1 TeV for m ψ /m S − 1 ∼ O(1). Finally, we have also estimated the limits from indirect detection experiments using gamma-rays and antiprotons as messengers.
We have found that the most stringent limits on the model stem from the non-observation of a gamma-ray flux correlated to the direction of dSphs. These limits are complementary to those from direct detection and collider experiments and exclude the region m ψ /m S 1.4, m S 150 GeV.