Cornering pseudoscalar-mediated dark matter with the LHC and cosmology

Models in which dark matter particles communicate with the visible sector through a pseudoscalar mediator are well-motivated both from a theoretical and from a phenomenological standpoint. With direct detection bounds being typically subleading in such scenarios, the main constraints stem either from collider searches for dark matter, or from indirect detection experiments. However, LHC searches for the mediator particles themselves can not only compete with -- or even supersede -- the reach of direct collider dark matter probes, but they can also test scenarios in which traditional monojet searches become irrelevant, especially when the mediator cannot decay on-shell into dark matter particles or its decay is suppressed. In this work we perform a detailed analysis of a pseudoscalar-mediated dark matter simplified model, taking into account a large set of collider constraints and concentrating on the parameter space regions favoured by cosmological and astrophysical data. We find that mediator masses above 100-200~GeV are essentially excluded by LHC searches in the case of large couplings to the top quark, while forthcoming collider and astrophysical measurements will further constrain the available parameter space.


Introduction
Searches for dark matter (DM) particles constitute one of the main physics objectives of the LHC, their prime signature being associated with the presence of missing transverse energy ( / E T ) in the collision final state. Traditionally conducted within specific dark matter models like supersymmetry that can give rise to a variety of final states involving / E T , LHC analyses have lately shifted towards more model-independent approaches [1]. One such approach is based on effective field theories (EFT) [2,3] which, despite its simplicity, can draw a misleading picture at LHC energies when TeV-scale new degrees of freedom cannot be integrated out [4][5][6]. Another approach relies on so-called simplified models, i.e. simple frameworks which extend the Standard Model (SM) by two particles, a dark matter candidate as well as a state that mediates the dark matter interactions with the visible sector [7][8][9]. Such frameworks palliate some of the drawbacks of the EFT approach, at the price of introducing only a handful of additional free parameters. Moreover, simplified models capture, with a minimal set of assumptions, some important features of more ultraviolet-complete (UV) theories and, perhaps more importantly, they can provide a (semi-)consistent framework in order to analyse the experimental results [10]. These two approaches are indeed complementary, since they explore different regions of the mediating particle's mass scale 1 and consistently set EFT limits can in principle be reinterpreted in any specific underlying model [12]. Independently of the theoretical framework, the most well studied model-independent dark matter signatures at the LHC have been the mono-X ones [13][14][15][16][17][18][19], in which a pair of dark matter particles is produced in association with a single energetic visible object. Amongst all mono-X searches, the monojet one has garnered the most attention given the relative magnitude of the strong coupling with respect to the electroweak one. This channel was shown to provide powerful constraints on dark matter models, especially in cases where direct detection experiments become inefficient [20,21]. Barring somehow singular kinematic configurations that can occur in freeze-out scenarios [22], the two most notable situations in which direct detection constraints fall short concern models in which the dark matter particles can be lighter than a few GeV, where direct detection searches suffer from recoil energy threshold limitations [23], or models in which the spin-independent dark matter-nucleon scattering cross section is suppressed due to the Lorentz structure of the underlying theory. This suppression can occur in models of axial-vector mediated fermionic dark matter or, which is the topic of this work, in scenarios featuring a pseudoscalar-mediated fermionic dark matter candidate in which the dark matter-nucleon scattering cross section is suppressed by the momentum transfer involved in the reaction.
If the pseudoscalar mediator, moreover, couples to the Standard Model fermions proportionally to their mass (an assumption motivated both by the success of Minimal Flavour Violation [24,25] and by ultraviolet considerations), it should be mostly produced through gluon fusion, analogously to Standard Model Higgs boson production. Gluon fusion processes typically exhibit an enriched jet activity compared to quark-antiquark annihilation ones. As a consequence, searches for final state signatures comprised of a multijet system accompanied by / E T could potentially lead to more stringent constraints than those originating from pure monojet searches, as first argued in Ref. [26] on the basis of 8 TeV LHC data. For this reason, more recent monojet LHC analyses at 13 TeV allow for an additional hadronic activity in their selection. In this work we investigate the constraints on the parameter space of a simplified pseudoscalar-mediated fermionic dark matter model stemming from both monojet and multijet plus / E T searches at the LHC. We quantify the impact of higher-order QCD corrections and we further study the bounds originating from complementary LHC dark matter searches such as the associated production of an invisibly-decaying pseudoscalar mediator with top or bottom quark pairs, also highlighting the future prospects for these searches.
Besides collider searches for dark matter in channels with large / E T , resonance searches by means of signatures made up solely of visible objects could be useful to constrain the properties of the mediator connecting the dark and visible sectors [27][28][29]. We therefore revisit LHC studies probing a potential new physics resonance decaying into a pair of tau leptons [30,31], photons [32] or top quarks [33]. Moreover, as we investigate scenarios with a large coupling of the mediator to top quarks, we also assess the sensitivity of the direct measurement of the top quark pair production cross section to the new physics parameter space. This observable turns out to play a key role even when a top-antitop pair is produced via an off-shell pseudoscalar exchange.
As the possible observation of enhanced missing energy signals (or, even more so, of a new resonance) does not guarantee their cosmological relevance, we compare the LHC constraints on the parameter space with the corresponding regions that are phenomenologically viable from a cosmological and astrophysical standpoint. Concretely, we investigate bounds arising from the dark matter relic density as well as from indirect searches such as the Fermi-LAT searches for dark matter-induced gamma-rays in Dwarf Spheroidal Galaxies [34] and for spectral features at the Galactic Centre [35], but also the AMS-02 searches for antiprotons [36]. Future prospects for indirect detection experiments are also discussed.
The paper is organised as follows: In section 2 we present the theoretical framework for our study. In section 3 we describe the various collider constraints that we consider, whereas section 4 is dedicated to dark matter astrophysical observables. Our results are presented in section 5 while in section 6 we provide our conclusions. Two appendices follow, discussing some technical aspects of our analysis.

Model description
We consider a new physics scenario in which the Standard Model field content is extended by a Majorana fermionic field χ of mass m χ , which plays the role of a dark matter candidate, and a pseudoscalar field A of mass m A , which mediates the interactions of χ with the Standard Model. Both particles are taken to be singlets under the Standard Model gauge group SU (3) c × SU (2) L × U (1) Y . Under these assumptions, and ignoring any cubic and quartic self-interaction of A, the part of the Lagrangian involving only the χ and A fields can be written down as where y χ denotes the strength of the interaction of the mediator with dark matter. Being a singlet under the Standard Model gauge group, A cannot couple to quarks and leptons through renormalisable gauge-invariant interactions. Hence, in order to parameterise the coupling of the dark sector with the visible one, we introduce the effective Lagrangian where the sums run over all up-type (f u ) and down-type (f d ) fermions respectively. The c u and c d coefficients parametrise the strength of the interactions between A and the SM fermions. In the spirit of Minimal Flavour Violation, we moreover take the corresponding operators to be proportional to the ratio of the SM fermion masses m f u,d over the vacuum expectation value of the Higgs field v. The Lagrangians of Eq. (2.1) and Eq. (2.2) will serve as a basis for the analysis that follows. In the above model description we have omitted, for the sake of simplicity, a quartic term involving bilinears of the A field and the SM Higgs doublet. Such a term could have interesting phenomenological consequences for dark matter, collider and Higgs phenomenology [37], and also on first order electroweak baryogenesis [38]. However, this falls beyond the scope of the present work.
It is also worth briefly commenting upon the potential UV origins of the Lagrangians (2.1) and (2.2) in order to motivate some of the parameter choices we will be adopting later on and to set the stage for the discussion that follows. The most straightforward UV completion of our setup would be in the framework of the two-Higgs doublet model (2HDM) or models involving even more extended scalar sectors. For example, in the context of type-II 2HDM, and denoting as usual by tan β = v 2 /v 1 the ratio of the two vacuum expectation values of the neutral CP-even components of the two Higgs doublets, we would obtain c u = cot β and c d = tan β. However, in such a scenario it is the Lagrangian (2.1) that would not be gauge-invariant. One solution could be to further introduce an additional scalar gauge singlet which mixes with the 2HDM Higgs doublets (an approach followed, e.g., in Ref. [39][40][41][42]) or, alternatively, to consider a pure 2HDM but extend the dark matter sector to a "bino-higgsino" or "bino-wino"-like fermion system [43,44].
In all these cases, additional interactions (e.g. with extra scalars) would arise at tree-level. Such interactions can introduce important phenomenological features which cannot be captured by the simple Lagrangians of Eq. (2.1) and Eq. (2.2). However, depending on the specificities of each potential generalisation of our framework, these additional features can be extremely model-dependent, rendering generic statements extremely hard (if possible) to extract. In this spirit, we adopt the simpler description provided by the Lagrangians of Eq. (2.1) and Eq. (2.2), cautioning the reader about subtleties that can appear in UV embeddings of our setup (see section 5.2). Similarly, we ignore constraints that can potentially arise from pre-cision electroweak tests or flavour observables as any realistic assessment of their impact would depend heavily upon the details of the UV embedding of the model. This has been addressed, for instance, in the framework of the 2HDM, in the work of Ref. [45]. Corrections to the S, T and U oblique parameters [46] or flavour constraints are nonetheless expected to be subleading for the A mass ranges under consideration [47]. In this work, we thus focus solely on direct probes of the dark matter particles χ and of the mediator A, whilst keeping track of the limitations of our simplified framework.
From a low-energy standpoint, the couplings y χ , c u and c d in Eq. (2.1) and Eq. (2.2) can take any numerical value (much like the new physics masses m A and m χ ) as long as perturbativity is respected. Throughout our study, results are shown for different discrete combinations of three out of these five parameters, the two others being varied freely. Our choices are mostly driven by phenomenological considerations while keeping in mind some model-building issues. In particular, we consider two distinct scenarios for the relative size of the c u and c d coefficients, namely which we refer to as "top-dominated scenarios" and that is referred to as a "bottom-dominated scenario". Interpreted, for example, in terms of a type-II 2HDM-like setup, the former case would correspond to tan β = 1 with standard (c u,d = 1) or somewhat enhanced (c u,d = 2) Yukawa couplings, whereas the latter to tan β = 10, again assuming slightly enhanced Yukawa couplings. However, these values for c u and c d have mostly been chosen because, as shown below, they allow for the exploration of different facets of the LHC phenomenology associated with our model. On the other hand, the dark matter relic abundance depends straightforwardly on the mediator mass m A . The latter is thus generically fixed when the dark matter phenomenology of the model is addressed. Conversely, as the LHC phenomenology of the model does not depend drastically on the dark matter mass itself (up to phase space considerations), but rather on its relation to the mediator mass, the discussion on the LHC constraints applicable to the model is performed within two setups. Either all coupling values are fixed and we vary the two new physics masses independently, or we fix the dark matter mass m χ and we vary its coupling to the mediator y χ along with m A .

LHC phenomenology
In order to check the compatibility of our new physics scenario with current LHC searches, we implement the model described in section 2 in the FeynRules pack-age [48] and export it under the UFO format [49] in order to make use of the MadGraph5 aMC@NLO platform [50] for the simulation of hard-scattering LHC collisions at centre-of-mass energies of 8, 13 and 14 TeV. When needed, these fixedorder results are matched with parton showers within the Pythia 6 environment [51] that is also employed for describing the hadronisation process, and we simulate the response of the LHC detectors with the Delphes 3 software [52]. When comparing our results with available experimental limits, we fold the cross sections computed through MadGraph5 aMC@NLO with appropriate K-factors to take into account non-simulated higher-order QCD corrections. For processes in which the leading-order contribution to the scattering amplitude arises at tree-level, as for the associated production of a top or bottom pair with dark matter, the corresponding K-factors are computed directly using our model implementation within the MadGraph5 aMC@NLO framework. In contrast, for processes whose lowest order diagrams are already at the one-loop level, the relevant K-factors are extracted from the literature, when available. These issues are further elaborated upon in the following paragraphs.
In order to systematise the discussion, we divide the presentation of the various LHC constraints between those involving invisible decays of the pseudoscalar mediator and those where the mediator decays into visible Standard Model objects.

Invisible mediator decay
Monojet and multijet plus missing transverse energy searches Monojet analyses are one of the primary probes for dark matter at the LHC, the targeted signature being characterised by the presence of a hard QCD jet recoiling almost back-to-back against a large amount of missing transverse momentum. Although they were originally designed to veto events in which any additional hadronic activity was present, it has been recently suggested that allowing for the presence of extra jets could improve the sensitivity of these searches, especially in the case where the partonic reaction is initiated by gluon fusion [26]. For this reason, both the ATLAS and CMS collaborations now include, in their monojet searches, signal regions populated by events involving more than one hard jet [20,21].
In this work we assess the compatibility of our scenario against the ATLAS monojet search with 3.2 fb −1 of integrated luminosity of proton-proton collisions at a centre-of-mass energy of 13 TeV [20]. We use a recasted version of this analysis [53] implemented in the MadAnalysis 5 framework [54][55][56] and available from the MadAnalysis Public Analysis Database 2 , following the procedure described in Appendix A for signal simulation. This analysis targets monojet events and contains various signal regions characterised by the considered amount of missing energy. Each region is associated with a different / E T selection threshold, the hardest se-lection corresponding to / E T > 700 GeV. In order to quantify the reach of such a monojet search for higher integrated luminosities, which potentially opens the door to more aggressive / E T thresholds, we define three additional signal regions in which the missing energy is required to be larger than 800, 1000 and 1200 GeV respectively. We extract the corresponding Standard Model background expectation from the official ATLAS / E T distributions that cover a missing transverse energy range extending up to 1200 GeV [20], and define the related uncertainty ∆N bkg on the number of expected background events N bkg as [57] ∆N 2 bkg = κ 1 N bkg 2 + κ 2 N bkg 2 with κ 1 = 1.5 and κ 2 = 0.043 .
In this expression, the first term represents the statistical error and the second the systematic one, the adopted values allowing us to adequately parameterise the AT-LAS results [20]. We apply this analysis on events for which the underlying matrix elements are allowed to contain up to one extra jet and that are merged, after parton showering, following the 'shower-k T ' merging scheme [58]. The first jet, already present at the level of the lowest jet multiplicity matrix element, is required to have a transverse momentum p T greater than 150 GeV to facilitate the accumulation of a higher statistics in the analysed signal regions, given the lowest experimental cut on the leading jet transverse momentum of 250 GeV. Our study shows that limits obtained with 300 and 3000 fb −1 of projected integrated luminosity are actually comparable with those obtained through the recasted version of the existing ATLAS analysis of 3.2 fb −1 of data, this lack of improvement being due to the ∼ 4% systematic uncertainty assumed in the background determination, a number which is unlikely to improve in the future. For this reason, we refrain ourselves from reinterpreting limits that could originate from more recent LHC analyses, such as the CMS monojet search with 12.9 fb −1 of integrated luminosity [21] which adopts a similar selection strategy and similar systematic uncertainty estimates as in our projections.
Motivated by the fact that for a gluon fusion process a higher jet multiplicity is expected in the final state, we have also examined the constraints that could arise from a supersymmetry-inspired multijet ATLAS search [59], basing our study on a recasted version of this analysis in the MadAnalysis 5 framework [60]. The fundamental difference between the monojet and multijet searches is that the latter involves a harder selection on the additional jets. Eventually, it turns out that the reduction in signal statistics outweighs the benefits of a more efficient background rejection, thus leading to slightly weaker limits. Going a step further, we have checked whether modifying a few selection cuts on the additional jets could improve upon the sensitivity of the monojet analysis by means of the multivariate analysis described in Appendix B. We have been unable to find any such improved set of cuts, thus concluding that under their present form, monojet searches appear to be optimal for targeting gluon-fusion-induced dark matter production processes. As limits from multijet searches turn out to be subleading with respect to the monojet-inspired ones, they are omitted in what follows.
ttA and bbA searches, with A → χχ The associated production of the pseudoscalar A with a pair of top or bottom quarks, i.e. the same topology as for the production of the Standard Model Higgs boson with a pair of third generation quarks, followed by the invisible decay A → χχ, could be an efficient probe of our model. This is especially true for top-dominated scenarios when the branching ratio BR(A → χχ) is substantial. Both the ATLAS and CMS collaborations have performed searches for an invisibly decaying spin-0 mediator particle produced in association with either a top or a bottom quark pair [61,62] 3 , so that these results can be reinterpreted as constraints on the scenarios studied in this work.
We consider results from the ATLAS search for an invisibly decaying heavy pseudoscalar mediator produced in association with a top quark pair in the singlelepton plus jets plus / E T channel [61]. This search analyses 13.2 fb −1 of proton-proton collisions at a centre-of-mass energy of 13 TeV and exclusion bounds at the 95% confidence level (CL) are presented as contours in the (m A , m χ ) plane for a fixed and common choice of the mediator coupling to the dark matter particle and to the top quark, c u = c d = y χ = 3.5. The experimental publication moreover includes results for smaller values of this common coupling setting in a restricted set of mediator and dark matter masses. Although the number of points is too limited to draw an exclusion contour, it is sufficient to check that this process does not constrain any further the parameter space of our model once other processes are accounted for, as discussed in section 5.1.
To deduce limits from the ATLAS search for all relevant masses, we compute, within the MadGraph5 aMC@NLO framework, the ttA associated production cross section at the next-to-leading order (NLO) in QCD, and further include leadingorder (LO) branching ratios for the decays of all particles. Furthermore, we consider the prospects of 300 fb −1 of LHC collisions at a centre-of-mass energy of 14 TeV. To this end, we import and reinterpret the projected upper limits on the signal strength, computed with respect to the case m χ = 1 GeV and c u = c d = y χ = 1, from Ref. [64]. Among the different investigations performed in this last study, we choose the shape-based results that assume a 20% uncertainty on the background estimation.
Conversely, in bottom-dominated scenarios, it is instead the associated production of the pseudoscalar mediator with a pair of bottom quarks that could potentially provide the strongest constraints. The ATLAS and CMS collaborations have performed analyses targeting this process for an integrated luminosity of 13.3 fb −1 [65] and 2.17 fb −1 [66] of LHC collisions at a centre-of-mass energy of 13 TeV, respectively. We make use of their findings to compare, as for the ttA case, the predicted signal cross sections with the excluded ones and derive exclusion bounds on the parameter space of our model.

Visible mediator decay
In addition to topologies where the mediator decays invisibly into a pair of dark matter particles, those involving decays into visible Standard Model states can be exploited to constrain our simplified model.
An important channel explored at the LHC is the production, either via gluon fusion or in association with a pair of bottom quarks, of a spin-0 mediator decaying into a pair of τ leptons. The CMS collaboration has performed an analysis targeting these topologies with an integrated luminosity of 12.9 fb −1 of LHC collisions at a centre-of-mass energy of 13 TeV. Final states featuring either 0, 1 or 2 hadronicallydecaying tau leptons have been equally considered [31], and the results have been presented as 95% CL upper limits on the production cross section of a heavy Higgs boson in the context of the Minimal Supersymmetric Standard Model. Both the gluon-fusion production channel and the associated bbA production mode have been constrained, and we confront these results to the predictions of our model. As no relevant information is provided explicitly, we assume that interference effects with the Standard Model have not been considered, and that the experimental results obtained for the case of a scalar mediator also hold in the pseudoscalar case.
For the gluon fusion channel, we multiply the LO cross sections returned by MadGraph5 aMC@NLO by N 3 LO A + N 3 LL K-factors for pseudoscalar mediators [67,68], so that the total rate includes the matching of approximate next-tonext-to-next-to-leading order predictions with the resummation of soft and collinear gluon radiation close to threshold at the next-to-next-to-next-to-leading logarithmic accuracy. In the case of the associated production of the mediator with a bottom quark pair, we consider instead NLO production rates multiplied by LO branching ratios.
Although the production of a Higgs boson in association with zero, one or two b-quarks and subsequently decaying into a pair of b-quarks or τ -leptons has also been explored at the Tevatron, the related sensitivity does not allow to probe pseudoscalar mediators lighter than 90 GeV [69,70]. No additional constraints can hence be deduced with respect to the ones extracted from the LHC results.

tt searches
One of the strongest constraints that can be imposed on our simplified model, especially in the case of top-dominated scenarios, comes from the study of new physics effects that can potentially appear within the production of a pair of top quarks. This corresponds to the pp → A ( * ) → tt process where the pseudoscalar can be either on-shell or off-shell. While dedicated LHC searches for tt resonances deeply probe the on-shell regime, the accurately measured tt production cross sections allow us to constrain the new physics parameters in the case of off-shell production as well.
For the latter possibility we consider several studies in which the total tt production cross section has been measured at the LHC, both for centre-of-mass energies of 8 and 13 TeV. The most precise 13 TeV measurement originates from the CMS analysis in the single-lepton plus jets channel, which yields a tt production cross section of [71] σ(tt) 13  on the size of the new physics contribution to the tt total production cross section. Similarly, the 8 TeV tt production cross section has been measured in several channels by the ATLAS [74][75][76] and CMS [77,78] collaborations. In order to extract constraints, we have started from the ATLAS and CMS measurements in the dileptonic mode [74,77], σ(tt) 8 = 242.9 ± 8.8 pb and σ(tt) 8 = 239 ± 13 pb , (3.5) that result in the same bound. We combine those numbers with the Standard Model expectation for m t = 172.5 GeV [73,79], which finally enables us to extract an upper bound at the 95% CL on the allowed size for the new physics contributions, all uncertainties being once again added in quadrature.
The 95% CL upper bounds are then compared to predictions within our new physics model. Interferences of the pp → A ( * ) → tt contribution with the Standard Model one should however be accounted for as they can be potentially large. Technically, they cannot be computed directly by MadGraph5 aMC@NLO and a trick must be employed. We first calculate the pp → A ( * ) → tt cross section within our model and next match it to the one obtained in the context of a dummy model where the heavy quark loop is approximated by a higher-dimensional effective operator involving the mediator and gluons. Within this second model, we then evaluate the interference with the SM tt background at LO, and we multiply this result by a K-factor that we take as the geometrical mean of the SM and new physics K-factors that are known separately. We find that the interference effects are important. They rise with increasing m A , are maximal near the tt threshold, and decrease for larger mediator masses. In particular, the interferences are often considerably larger than the new physics contribution itself when the mediator mass lies right below the threshold region. They must thus be included before evaluating the constraints on the model.
We furthermore evaluate the potential constraints arising from resonance tt searches. We use the Run I ATLAS results obtained for an integrated luminosity of 20.3 fb −1 of 8 TeV collisions [33], assuming that the quoted limits for a scalar mediator hold in our setup.

Diphoton searches
The pseudoscalar mediator A couples to the ordinary fermions in a similar manner as the Standard Model Higgs boson, proportionally to their mass. This, in turn, implies that it can decay into a pair of photons through a loop involving top or bottom quarks, with the former contribution typically dominating over the latter unless c d c u as in our bottom-dominated scenario. The fact that the Standard Model Higgs boson was first observed in the γγ channel motivates us to study constraints stemming from searches for diphoton resonances. The partial decay width of a pseudoscalar A into a photon pair is given by [80] where G f and α are the Fermi and fine-structure constants respectively, N c is the number of colours associated with each of the fermions running in the loop, Q f their electric charge and the sum runs over all Standard Model fermion species. Moreover, τ f = m 2 A /4m 2 f and the loop function A A 1/2 is given by (3.10) We employ the limits presented by the ATLAS collaboration in Ref.
[32] where 15.4 fb −1 of 13 TeV collision data is used. Based on the diphoton invariant mass distribution, this model-independent analysis targets spin-0 resonances with masses as low as about 200 GeV, and the results are presented as upper limits on the production cross section times branching ratio for different choices of the resonance width. Similarly, we also make use of the corresponding 8 TeV ATLAS results [81] that extend the covered mass range down to about 65 GeV.
For the entire considered parameter range, we have verified that the total decay width of the mediator is always small enough so that the narrow width approximation holds. This allows us to directly use the ATLAS limits and compare them with our predictions at the N 3 LO A + N 3 LL accuracy for what concerns the mediator production rate, the corresponding branching ratio into a photon pair being evaluated at LO. The interference effects with the continuum SM background have again been ignored here.

Other constraints
Several other collider searches could in principle also constrain the class of dark matter scenarios under consideration. Amongst the visible mediator decay processes, we have checked that mediator-induced dijet production does not enforce any relevant constraint on the parameter space [82], as do other mono-X searches such as the monophoton ones [83,84]. Another possibility are constraints arising from searches for the presence of pseudoscalars in the decays of the SM Higgs boson, which could be relevant for small enough mediator masses. Such searches have been performed, e.g., in the 4τ channel [85,86]. However, given the restricted form of the scalar potential (see section 2), these limits do not apply in our case.
As regarding constraints arising from LEP data, searches were conducted for the associated production of the mediator along with a pair of bottom quarks, with the subsequent A → bb or A → τ + τ − decay, and were used to derive constraints on pseudoscalars for masses lying in the [5,50] GeV window. The obtained limits, at the 95% CL, impose that for m A = 5 (50) GeV in the τ + τ − final state and for m A = 12 (50) GeV in the bb final state [87]. Due to its larger branching ratio, the bb channel leads to the most stringent upper limit on c d , but the derived constraints are evaded as we consider larger mediator masses not reachable at LEP. Constraints from precision measurements on simplified models should be interpreted with care as additional contributions or cancellations can occur within a UV-complete theory. We therefore only comment on their potential impact. Pseudoscalar contributions to the gauge boson self-energies appear at the two-loop level and we therefore only expect weak constraints on the model from precision electroweak measurements. Constrains on light pseudoscalars (m A < 10 GeV) can be obtained from flavour physics, in particular from measurements of the B s → µ + µ − decay. Imposing the loose requirement that the pseudoscalar contribution does not exceed the SM expectation leads to the bound [47] √ (3.14) The full range of parameters considered in this work is thus not affected.

Dark Matter phenomenology
The observation of missing energy signals or the detection of a new resonance at the LHC only provide little information concerning the cosmological relevance of the underlying physics. In this spirit, we wish to confront the constraints arising from LHC searches with those stemming from DM-related experiments, enforcing in this way that the properties of the χ particle do not challenge the cosmological and astrophysical observations. With large enough couplings such as to yield observable rates at the LHC, our dark matter candidate thermalises in the early universe and its present abundance can be computed in the framework of a standard thermal freeze-out. As the LHC phenomenology connected to CP -even and CP -odd scalar mediators is largely identical, the dark matter observables are actually the ones that have dictated our initial choice of considering pseudoscalar rather than scalar mediators. CP -even mediators are indeed essentially excluded by the results of direct detection experiments such as LUX [23] once they are combined with relic density requirements, under the condition that caveats like the invocation of mechanisms such as entropy injection that dilutes the DM abundance are ignored [90]. In the following subsections, we proceed to a brief description of the experimental constraints used in this work. All dark matter observables hereafter have been computed with the micrOMEGAs code [91][92][93] that relies on CalcHEP [94] model files obtained from the implementation of the model in FeynRules [95].

Relic abundance
The abundance of dark matter in the universe today has been precisely measured by the Planck mission [96]. At present, its central value reads

Indirect detection
With direct detection constraints being largely inefficient to probe the parameter space of our dark mater scenario, the main astrophysical constraints originate from DM indirect detection and, in particular, from measurements of the continuum gamma-ray spectrum from dwarf spheroidal galaxies, from searches for spectral features at the Galactic Centre and from measurements of cosmic ray antiproton fluxes.
Dwarf spheroidal galaxies (dSphs) are dark matter-dominated objects, a property which makes them very clean targets to test dark matter interactions with the visible sector. The most recent search for gamma-rays in dSphs has been performed by the Fermi-LAT collaboration [34], and experimental bounds on the thermally averaged DM annihilation cross section have been derived under the assumption of a 100% annihilation into given Standard Model channels. However, in our model, dark matter can annihilate into any SM fermion pair, if kinematically allowed, with a specific weight ω that depends on the choice for the couplings c u and c d in Eq. (2.2). In order to account for the presence of several annihilation modes contributing to the gamma-ray signal, we recast the experimental bound on the thermally-averaged total annihilation cross section as where we sum upon all possible annihilation modes, where c i generically stands for c u and c d , and σv exp j are the reported experimental 90% CL limits on the annihilation cross section assuming a 100% annihilation into the channel j. The limit deduced from Eq. (4.4) is then compared with the σv value predicted by the model. As DM annihilation into a pair of pseudoscalars is p-wave-suppressed, as mentioned above, this channel is always subdominant and thus safely ignored from our analysis of the DM indirect detection bounds. In a similar fashion, we moreover report the expected 15-year exclusion bounds that could be obtained assuming data issued from 60 dSphs, following the projections of the Fermi-LAT collaboration [97]. We, however, do not attempt to fit the Galactic Center excess observed by FermiLAT [98,99], which has been considered in the context of pseudoscalar mediated simplified models in Refs. [100][101][102] and EFT [103].
On different lines, cosmic ray antiprotons constitute at present one of the most important channels for indirect dark matter searches. Recently, the AMS-02 experiment has released its results on the cosmic ray antiproton-to-proton ratio with an improved statistical precision [36], which has recently allowed to derive an approximate 2σ limit on the thermally-averaged DM annihilation cross section in the bb channel [104]. We assume that this last limit is representative for all quark flavours [105] and then proceed to extract the maximum allowed annihilation cross section according to the expression of Eq. (4.4). The ensuing constraints however depend strongly on astrophysical assumptions such as the cosmic ray propagation model and the DM density profile. We use as a benchmark an Einasto DM halo profile [106] and consider the so-called MED propagation model. As we will show in the following, for this benchmark the antiproton channel exclusion power is largely subdominant with respect to the sensitivity of the dSphs data. Given that dSphs constraints are generically considered to be more robust we will, thus, refrain from showing antiproton bounds adopting different astrophysical assumptions. We note, however, that as shown in [104], opting for a MAX propagation model would amount to bounds on σv which would be stronger by roughly one order of magnitude. A similar conclusion can be drawn from the recent analysis performed in [107], with the exception of a dark matter mass range around 100 GeV where, depending on the whether or not the low-energy part of the AMS-02 antiproton data is taken into account, this group reports in general weaker limits. Besides, varying the assumed dark matter halo profile is known to only mildly affect the relevant constraints [104].
Finally, the coupling of the mediator A to the Standard Model quarks can induce, at the one-loop level, DM annihilation into a pair of photons that could be observed as gamma-ray lines. The most stringent bounds on such a mechanism are issued from the Fermi-LAT searches for spectral features at the Galactic Centre [35]. We estimate the related constraints on our model by first matching the expression of Eq. (3.8) for the mediator diphoton decay width to the one obtained when relying on an effective interaction Lagrangian, where F andF respectively stand for the photon field strength tensor and its dual.
The m A dependence of the form factor A A 1/2 in Eq. (3.8) is however replaced by a 2m χ dependence, the latter being the relevant energy scale of the process instead of the mediator mass. The DM annihilation cross section into a photon pair is then computed in the resulting effective field theory, with the scale Λ γ being appropriately fixed by the above procedure.

Other constraints
Similarly to indirect detection, Cosmic Microwave Background data can give rise to a priori relevant constraints as well. Even if DM does not significantly annihilate directly into electrons in our model, other highly energetic annihilation products can heat and ionise the intergalactic medium, although with less efficiency, and thus affect the last scattering surface which is measured with high precision by the Planck collaboration. Taking the most updated limits from Ref. [108], we have verified that any related constraint is roughly one order of magnitude weaker than those detailed in the previous sections.
As already mentioned, the bounds on our model that can be obtained from DM direct detection data are very weak as the effective operator that describes the DM-nucleus interactions is momentum-suppressed. The current limits obtained by the Xenon collaboration [109] on several classes of operators (and in particular the so-called operator O 6 [110,111] that could be relevant in this work) do not further constrain the parameter space under consideration.

Dark-matter-favoured regions of the parameter space
Before closing this section and analysing the impact of LHC searches on our model, we first determine which regions of the parameter space are favoured by DM considerations. For given values of the mediator mass and of its couplings to the Standard Model fermions, accommodating the correct relic density leads to an m χ -dependent lower bound on the mediator coupling to dark matter y χ . The strength of this coupling has to be large for light DM, is minimum for m χ ≈ m A /2 and increases again for larger DM masses until a new efficient annihilation channel opens up at the topquark threshold m χ ∼ m t . This in turn requires a small y χ value to saturate the Planck bound. Moreover, although we mainly restrict ourselves to the parameter space region in which m χ < m A , a new annihilation channel χχ → AA opens once the m χ = m A threshold is crossed, which may dominate the total DM annihilation depending on the parameters of the model. In this case, the relic density predictions become independent, if the narrow-width approximation holds, of the coupling between the mediator A and the Standard Model.
In Figure 1, we present, as exclusion contours, all the constraints on our model that have been discussed in this section. We fix the mediator mass to m A = 250 GeV and coupling strengths to the SM to c u = c d = 2, and show results in the (m χ , y χ ) plane. The various combinations of parameters that can account for the entire dark matter abundance in the universe according to standard thermal freeze-out are represented by a black line. We can observe that the limits originating from the Fermi-LAT observations of dwarf spheroidal galaxies, shown as a solid red line, are the most powerful ones for DM masses below the electroweak scale. These constraints exclude the generic s-wave DM annihilation cross section favoured by Planck when DM can only annihilate into light quarks or tau leptons. Consequently, the Fermi-LAT data strongly disfavour the parameter space region compatible with the Planck lower bound for the relic density for m χ < 70 GeV. For larger DM masses, the constraints become in contrast weaker and there is always a set of couplings for which both Planck and dSphs constraints can be satisfied simultaneously, although this region can be very narrow. Similar results are found for other values of c u = c d as long as y χ × c f is kept constant. The exact location of the allowed region moreover depends on the choice of the mediator mass, but a region allowed by all DM constraints can always be found for pseudoscalar masses in the 10-1000 GeV range provided c u and c d are both not too suppressed. The projected 15-year Fermi-LAT limit, depicted by a red dashed line, demonstrates that a conclusive statement about the viability of the scenarios under consideration will be feasible. Almost the entire parameter space compatible with Planck data is expected to be excluded, except for a small resonant region around m χ ≈ m A /2, or if the DM particles are very heavy.
We also see from Figure 1  data for the entire considered mass range. They however supersede the Fermi-LAT bounds for large couplings y χ ≈ 2 and light DM with m χ < 50 GeV. These results are nonetheless not as robust as they include a substantial amount of assumptions, in particular concerning the cosmic-ray propagation model. Besides, constraints from gamma-ray line searches are always subleading when compared to gammaray continuum analyses, as the higher sensitivity does not compensate the loopsuppression factor entering the DM annihilation cross section into photon pairs in our model. Given these findings, we consequently omit, in the following, all constraints stemming from antiprotons and gamma-ray lines and focus instead on the limits and projections originating from searches for the gamma-ray continuum in the dSphs.

Results and discussion
Having presented the collider searches we consider in our analysis and sketched the parameter space region that is favoured by dark matter observables, we now estimate the impact of both the LHC and the astrophysical constraints on our model.

Fixed couplings
As a first step, following standard practice for the presentation of LHC dark matter searches [1,7,112], we fix the mediator couplings to the visible and dark sectors and study the interplay of collider and DM observables in the (m A , m χ ) plane. LHC probes are more effective when the pseudoscalar couples strongly to the Standard Model. In particular, we have found that the tt searches almost single-handedly exclude the case c u = c d = 3 for m A even below 100 GeV and up to more than 1 TeV. We therefore choose to show our results setting the couplings of the pseudoscalar to ordinary fermions to lower values, picking as a benchmark scenario whereas we fix the coupling to dark matter particles to Our results are shown in Figure 2, where we highlight the (m χ , m A ) combinations saturating the Planck bound (black solid lines), along with the parameter space regions excluded by the Fermi dSphs observations (red solid lines) and 15-year projections (red dashed line), monojet searches (hatched green region), A → τ + τ − searches (grey region), diphoton resonance searches (blue region) and tt searches at 8 TeV (hatched grey region). There are two regions preferred from a dark matter standpoint. The first one corresponds to m A > 300 GeV (the area above the black and below the red lines towards the right-hand side in Figure 2) and the other to lower values of m A with m χ ≈ 150 GeV (the area between the black and red lines towards the top-left of the same figure). For m A 100 GeV, dark matter is underabundant for small DM masses and the relic density increases with m χ until it reaches the Planck limit around m χ ∼ 110 GeV. Then, as soon as the tt annihilation channel opens up, the DM annihilation cross section increases significantly and χ becomes underabundant again. Besides, due to this sudden increase in σv , m χ values above m t are excluded from Fermi-LAT. However, because of the finite velocity of dark matter particles in the early Universe, this increase occurs at a somewhat lower value of m χ for the annihilation cross section during freeze-out (concretely, around m χ 155 GeV). This explains why dark matter masses around 160 GeV are more difficult to probe with Fermi-LAT, as can be seen in Figure 2: they correspond to situations where the observed dark matter abundance is obtained through annihilation into tt pairs in the early Universe, with this process being inefficient at present times. In any case, the 15-year projected limits of Fermi-LAT can fully exclude the parameter space regions of moderate pseudoscalar masses. Because of the strong coupling to top quarks, the region with a heavy pseudoscalar is severely constrained by the LHC. The 13 TeV tt total cross section measurement excludes the considered model configuration for pseudoscalar masses m A between 310 GeV and 410 GeV 4 , while the corresponding 8 TeV one further pushes the lower limit to m A = 280 GeV and the upper one to 430 GeV. The strongest constraint beyond m A = 400 GeV comes from the 8 TeV tt resonance search which excludes the considered benchmarks for mediator masses ranging up to ∼ 800 GeV.
The process Att where A decays invisibly does not provide additional constraints: we have found that only DM masses m χ < m A /2 are excluded at the 95% CL as long as m A < 200 GeV, while for m A > 300 GeV all parameter space points are allowed 5 .
The monojet and multijet constraints are only efficient for m A > 2m χ and their reach extends up to m A = 320 GeV, mostly probing Planck-preferred regions which are already challenged by Fermi-LAT. The searches for diphoton resonances also constrain the region 100 m A 350 GeV. The lower value is due on one hand to the suppression of the diphoton branching ratio for small m A values and on the other hand to the reduced sensitivity of the LHC to light diphoton resonances. The upper value is a consequence of the fact that once the tt threshold is reached, the decay width of A into top-quark pairs dominates and the diphoton channel constraints disappear. The non-trivial m A -dependence of these limits in the intermediate region is due to a combination of two effects: generically, once the decay A → χχ becomes kinematically allowed the branching ratio into photon pairs decreases substantially and the diphoton constraints vanish. However, as m A approaches 2m t the loop form-factor (3.9) entering the diphoton decay width (3.8) becomes maximal and the constraints become stronger, until the A → tt channel opens up. Besides, direct searches for pseudoscalar decays in the τ + τ − final state cover m A masses in the [120-320] GeV range, provided the branching ratio into a ditau system is not too suppressed by the invisible partial width. In this channel, extending the search to lower masses will be required in order to probe the dark matter-allowed region for light pseudoscalars. We remind that all collider limits presented here correspond to a 95% CL exclusion.
Reducing the coupling of the pseudoscalar to fermions (and in particular c u ) would considerably weaken the LHC constraints and especially the tt one as will be seen below. The limits from the diphoton channel also become weaker, since the loopinduced production cross section scales as c 2 u , with the bottom quark contribution to the process being subleading.

Fixed SM couplings and DM mass -benchmarks
Next, we investigate the impact of the various constraints on the parameter space for fixed values of the DM mass. We choose and, as discussed in section 2, three sets of couplings as in the bottom-dominated scenario. The results are projected in the (m A , y χ ) plane.

Fixed SM couplings and DM mass -scenario S1
We consider first the scenario S1 in which c u = c d = 2, for which our results are summarised in Figure 3. As mentioned before, tt searches place severe constraints on this scenario, in particular due to the measurement of the tt cross section at 8 TeV that entirely excludes the region m A > 280 GeV. We stress again that interference effects are important here, and fairly independent of the value of y χ . The remainder of the DM-favoured region for m A > 2m χ is also excluded by various searches. The gg → A → τ + τ − search eliminates the region from m A = 130 GeV up to 220 GeV for all values of y χ . When the pseudoscalar can decay invisibly, its branching ratio into τ pairs is suppressed, inducing a non-trivial dependence of the exclusion bounds on y χ . For y χ < 0.5, the limits extend to larger pseudoscalar masses reaching m A = 350 GeV for y χ = 0.3. They sharply drop for masses above the tt threshold, where no exclusion can be obtained from this channel. These limits can potentially be extended to lower values of m A , under the condition the experimental results become publicly available.
Diphoton probes cover the region 120 GeV < m A < 380 GeV, where the general behaviour of these constraints can be understood following a similar line of reasoning as before: the lower value is due to the suppression of the diphoton branching ratio for smaller m A and to the reduced LHC sensitivity. The upper value results from the same reduced production cross section combined with a reduced diphoton branching ratio given the large partial width into top pairs. Again, once m A becomes larger than 2m χ , a large value of y χ implies an increased total width and a reduced branching ratio into photons, with the exclusion limits featuring a non-trivial dependence both on y χ and on m A . More precisely, as soon as the decay A → χχ becomes kinematically allowed, we observe a sharp drop in the y χ values that are reachable due to the subsequent decrease in the diphoton branching ratio. On the other hand, as m A increases further, the constraints become stronger due to the maximisation of the loop form-factor of Eq. (3.9) entering the diphoton decay width of Eq.
The monojet search possesses a similar sensitivity on the mediator mass, the 220 GeV < m A < 340 GeV mass range being covered, but the results also depend on y χ as the search becomes ineffective when the pseudoscalar coupling to DM is too small. In this sense, it is complementary to the diphoton search.
The projection for 300 fb −1 of future LHC collisions at a centre-of-mass energy of 14 TeV shows that ttA probes (with A → χχ) will allow for the coverage of the parameter region currently probed by monojet-like searches for mediator masses ranging up to m A = 350 GeV and that the results depend on the value of y χ below the tt threshold. For smaller y χ values of, for instance, 0.2, the projected exclusion extends only to lower masses.
In summary, the only DM-favoured region that survives the various LHC bounds lies below m A = 130 GeV and conclusive statements on the viability of this scenario will be achievable in the future, thanks to the Fermi-LAT observations.

Fixed SM couplings and DM mass -scenario S2
For weaker couplings to the SM particles, the LHC constraints get significantly relaxed. The case c u = c d = 1 (scenario S2) is illustrated in Figure 4, where we for instance observe that searches with tt probes become less efficient. There is only a small corner of the S2 parameter space, in which 400 GeV < m A < 450 GeV, that is excluded by 8 TeV tt resonance searches. In addition, configurations for which  Finally, the monojet search in contrast barely excludes a very narrow mass range around m A = 200 GeV for coupling y χ 1.3, which is not shown in the Figure. In summary, most of the S2 parameter space region favoured by DM (between the black and red lines) is consistent with all current LHC constraints. The mass region 200 GeV < m A < 320 GeV is expected to be covered by the future results of the ttA LHC analyses of 14 TeV collisions, as long as y χ > 0.1. Fermi-LAT will, besides, be able to exclude essentially the entire parameter space with 15 years of data acquisition.

Fixed SM couplings and DM mass -scenario S3
We finally consider a case in which the coupling to b-quarks is enhanced and the one to t-quarks is suppressed (scenario S3). As for the other considered cases, we have found two viable DM regions, one for m A > 350 GeV and another (very narrow) one for m A < 90 GeV. Clearly, no limits can be derived from the tt or ttA processes. The main LHC constraints instead come from searches in the bbA with A → τ + τ − channel. Pseudoscalar masses between m A = 90 GeV and 450 GeV are fully excluded for y χ < 1.5. The exclusion bounds extend to larger pseudoscalar masses but for smaller dark matter-mediator couplings, reaching up to m A = 660 GeV for y χ 0.1. There is also a narrow excluded region between m A = 90 GeV and 100 GeV from  the pseudoscalar (produced in the gluon fusion mode) search in τ + τ − . Again, the absence of data below m A = 90 GeV makes it difficult to predict the viability of the considered scenarios in this region. Additional constraints could be obtained from the search for Abb with A decaying to a pair of dark matter particles. The limit extracted from Ref. [66] for m χ = 100 GeV is not available, however it should be weaker than the quoted limit for m χ = 1 GeV, for a given m A > 2m χ . We found that the regions excluded by this constraint would lie in the range 200 GeV< m A < 300 GeV for y χ 1.5 and are, thus, anyway covered by the bbA, A → τ + τ − searches.
In summary, the DM-favoured regions in the bottom-dominated scenario survive the LHC constraints for m A > 620 GeV and m A < 90 GeV. As for our other scenarios, we expect both these regions to be probed by Fermi-LAT in the future.

Discussion
In the previous subsection we presented results fixing three out of the five free parameters of our model at a time. In order to ascertain the generality of these results, we comment below on the behaviour of the limits for different choices of the dark matter mass and mediator couplings. On general grounds, LHC constraints are mostly important for large couplings of the pseudoscalar to top quarks, at least for the top-dominated scenarios S1 and S2. As illustrated in Figure 1, DM masses below ∼ 70 GeV are strongly constrained from a combination of Fermi-LAT and Planck results for fixed (and actually rather large) values of the couplings c u and c d . Moreover, decreasing the DM mass generically shifts the Planck-compatible region towards higher values of y χ pushing it quickly into the large coupling regime. As for LHC constraints, both the monojet and ttA exclusion bounds are expected to be extended to lower values of m A (down to m A ∼ 2m χ ). Alternatively, increasing the DM mass shifts the Planck-favoured region associated with m A > 2m χ towards smaller couplings y χ for a fixed pseudoscalar mass, and therefore this region is constrained mainly from the tt channel (see e.g. Figure 3). Other constraints are mostly limited to the region with m A < 2m t and are largely independent of the other parameters of the model -for example, the precise value of the DM mass will only shift the lower limit of the exclusion bounds from invisible channels, which is anyway confined to m A > 2m χ . Finally, the limits coming from the Fermi-LAT dSphs searches become weaker for heavier dark matter, hence, the DM-favoured region will become more extended as m χ increases.
The results for the bottom-dominated case show that the main LHC constraints do not involve invisible final states, and are therefore mostly insensitive to the DM mass. However, the limits originating from A → τ + τ − searches will become stronger when the A →χχ decay is not kinematically open, thus covering the low m A values which are missed by the bbA, A → τ + τ − searches. For example, we have checked that increasing the value of c d to 30 would amount to an exclusion of pseudoscalar masses up to m A ∼ 860 GeV.
As a final remark, we should stress again the limitations of simplified models. In particular, the relic density predictions can vary substantially in concrete UV completions of our model due to the presence of additional annihilation channels involving, for example, extra scalars 6 . Such additional contributions to DM annihilation in general affect the relic density and the gamma-ray flux in the same direction and, hence, the parameter space regions allowed by both Planck and Fermi-LAT are expected to remain narrow. The main exception to this rule is connected to co-annihilations, which only affect the relic density and therefore reduce the region where DM is overabundant. As for LHC limits, for m A < 2m t , the presence of additional decay channels can indeed substantially modify our results. As long as m A > 2m t , we however expect the limits obtained in this work to remain qualitatively valid (provided that the narrow-width approximation holds). Exceptions do exist, for instance in scenarios where additional contributions to the tt total cross section interfere with the SM and/or pseudoscalar ones. The situation becomes, of course, more involved once we consider the evolution of the interplay between cosmological/astrophysical constraints and collider ones, a case in which the hierarchy between all masses involved in potential UV completions of our simplified model should be taken into account.

Summary and outlook
From a combination of searches for particles decaying into Standard Model and/or invisible final states, we have shown that the LHC provides strong constraints on models of fermionic dark matter coupled to the visible sector through a pseudoscalar mediator A. The constraints apply mainly when at least one of the mediator couplings to quarks is of O(1) and cover a substantial fraction of the DM-favoured region of the parameter space.
When the mediator possesses a sizeable coupling to top quarks, LHC dark matter searches based on invisible channels are only relevant for m A 350 GeV as for higher masses its total width is dominated by the decay into tt. Direct searches for the mediator in channels involving top quark pairs overcome this restriction and, as we have shown, allow the LHC to probe pseudoscalar masses up to 800 GeV, depending on the exact value of the relevant couplings. A dedicated search for resonances in the tt channel with increased luminosity is, thus, expected to significantly extend the LHC reach to larger mediator masses and to smaller values of its couplings to quarks. Moreover, searches for an invisibly decaying pseudoscalar produced in association with top quarks will be able to probe a significant fraction of the DM-favoured region when the mediator mass is in the range 2m χ < m A < 2m t . Our findings further show that the projected limits obtained from the ttA channel are complementary to the bounds stemming from searches for diphoton resonances as well as to the corresponding bounds from traditional monojet searches. Besides, pseudoscalar masses below 100-200 GeV survive all constraints, where the exact value depends on the choice of the DM mass and of the mediator couplings. The reason is either a lack of data in the region of interest, e.g. in the ditau channel, or a suppressed branching ratio of the mediator, e.g. in the case of diphotons. Thus, most of the DM-favoured regions at mediator masses of O(100) GeV or lower remain allowed once LHC constraints are imposed. However, if the mediator couples strongly to down-type fermions (in the so-called bottom-dominated scenario, S3), this mass range could be probed at the future LHC runs through searches for light pseudoscalars in the bbA channel, with the mediator decaying into a pair of taus. The light pseudoscalar scenario can also easily be probed at a future electron-positron collider in either the τ + τ − or bb mode, as for such masses the associated cross section exceeds 20 fb [113].
On the side of indirect detection, searches for gamma-rays from Dwarf Spheroidal galaxies will allow Fermi-LAT to cover most of the DM-favoured region after 15 years of data acquisition. One known exception is the DM resonance region where m χ ≈ m A /2. In this case, the mediator couplings either to dark matter or to the Standard Model particles (or, eventually, to both) have to be very suppressed in order to saturate the relic density bound. In a narrow mass range, the total thermally averaged annihilation cross section σv at galactic velocities indeed turns out to be much suppressed compared to σv in the early universe, rendering even future indirect detection observations irrelevant. Although it could be argued that this parameter space region is fine-tuned, it generically occurs in all models where DM annihilation proceeds through an s-channel resonance. The LHC, however, can (at least partly) probe the DM resonance region through searches for the mediator decaying into visible final states such as top or tau pairs, provided its couplings with the Standard Model particles are not suppressed. This feature is yet another illustration of the complementarity between astroparticle and collider searches for dark matter. Satyaki Bhattacharya for help with multivariate analyses, to Narayan Rana for clarifications about K-factors adopted for gluon-fusion processes and to Nicolas Berger for ascertaining the fact that a combined Gaussian uncertainty is conservative even if the actual distributions deviate from exact Gaussian shapes. This work was supported in part by the French ANR project DMAstro-LHC (ANR-12-BS05-0006), by the Investissements d'avenir Labex ENIGMASS, by French state funds managed by the Agence Nationale de la Recherche (ANR) in the context of the LABEX ILP (ANR-11-IDEX-0004-02, ANR-10-LABX-63), and by the Research Executive Agency (REA) of the European Union under the Grant Agreement PITN-GA2012-316704 (HiggsTools). SB acknowledges support from the Indo-French LIA THEP (Theoretical High Energy Physics) of the CNRS.

A Merging and matching
In order to accurately simulate the monojet and multijet processes described in section 3.1, we have generated hard scattering events for the pp → A process with matrix elements containing up to one extra jet. These fixed-order results have been matched with the parton shower infrastructure of Pythia 6 and then combined following the 'shower-k T ' merging scheme [58]. The merging parameters Q cut (at the matrix-element level) and the merging scale Q have been fixed to the common m A -dependent value reported in Table 1.
We have moreover verified the stability of our results with respect to the case in which up to two extra jets are allowed at the matrix element level. Due do the prohibitive computational cost of performing this task with the full model described in section 2, we have opted for a simpler effective field theory model described by the Lagrangian where B µν , W µν and G µν are the U (1) Y , SU (2) L and SU (3) c field strength tensors respectively, andB µν ,W µν andG µν their duals. Moreover, g 1 , g 2 and g 3 are the hypercharge, weak and strong coupling constants. We have verified that the two procedures yield comparable results once the experimental selections for the monojet and multijet analyses are applied for a selected set of scenarios, justifying our choice of simulating hard scattering events with up to only one jet in the matrix element.

B MVA analysis
Having employed the ATLAS monojet-like analysis in order to constrain our model, we further attempted to check if any kind of improvement is possible in this regard.
To this end, we again relied on the EFT scenario described by the Lagrangian of Eq. (A.1). As discussed in section 3.1, the ATLAS monojet-like analysis at 13 TeV relies on a strategy that requires the presence of a hard jet with a transverse momentum p T > 250 GeV and of at most four additional jets with p T > 30 GeV. Seven inclusive and six exclusive signal regions have been studied in the ATLAS search, all defined by a multitude of selection cuts (some more details on these issues have been given in section 3.1). In order to study whether any further sensitivity can be gained, we performed a multivariate analysis (MVA) relying on a boosted decision tree (BDT). We made use of the TMVA framework [114] and considered 14 kinematic variables, namely the missing transverse energy / E T , the p T , the pseudorapidity and azimuthal angle of the two leading jets, the angular separation in azimuth between the two leading jets as well as between each of them and the missing momentum, the angular distance in the transverse plane between the two leading jets, the transverse mass reconstructed from each of the two leading jets and the missing transverse momentum, the angular separation in azimuth between the two leading jets as well as between each of them and the missing momentum, the angular distance in the transverse plane between the two leading jets, the transverse mass reconstructed from each of the two leading jets and the missing transverse momentum and the aplanarity [115].
Throughout our MVA analysis, we carefully treated the issue of the overtraining and we internally performed the Kolmogorov-Smirnov (KS) test. Whereas in principle, the KS probability associated with both the signal and the background needs to lie between 0.1 and 0.9, a critical KS probability larger than 0.01 [116] and not exhibiting oscillations is also acceptable.
For our BDT analysis, we have only considered the two dominant backgrounds consisting of invisible Z-boson plus jets and leptonic W -boson plus jets events, carefully taking into account their individual weights. For our first signal sample, we imposed a very loose selection cut before using the BDT, constraining the p T of the leading jet to be larger than 150 GeV and the number of jets to be at least 2. Choosing a benchmark setup with Λ 3 = 700 GeV, m A = 780 GeV and m χ = 190 GeV, we obtain a number of signal and background events of N S = 645 and N B = 9124 respectively for an integrated luminosity of 3.2 fb −1 . The BDT analysis yields a significance σ 0 = N S / √ N S + N B = 6.49 for a ratio N S /N B of about 7%, assuming zero systematic uncertainties. Upon considering a flat 10% systematic uncertainty on the SM backgrounds, the significance drops to σ 10 = 0.70, which is almost exactly the same as the one obtained from the best signal region in the cut-based analysis of ATLAS. As a final check, we relaxed the initial p T selection cut on the first jet to 80 GeV and observed that the N S /N B ratio drops to about 3% and σ 0 to 6.17.
These findings lead us to the conclusion that the ATLAS cut-based analysis is already highly optimised and, as our multivariate analysis has shown, further improving upon it appears to be highly non-trivial.