A comprehensive approach to dark matter studies: exploration of simplified top-philic models

Studies of dark matter lie at the interface of collider physics, astrophysics and cosmology. Constraining models featuring dark matter candidates entails the capability to provide accurate predictions for large sets of observables and compare them to a wide spectrum of data. We present a framework which, starting from a model lagrangian, allows one to consistently and systematically make predictions, as well as to confront those predictions with a multitude of experimental results. As an application, we consider a class of simplified dark matter models where a scalar mediator couples only to the top quark and a fermionic dark sector (i.e. the simplified top-philic dark matter model). We study in detail the complementarity of relic density, direct/indirect detection and collider searches in constraining the multi-dimensional model parameter space, and efficiently identify regions where individual approaches to dark matter detection provide the most stringent bounds. In the context of collider studies of dark matter, we point out the complementarity of LHC searches in probing different regions of the model parameter space with final states involving top quarks, photons, jets and/or missing energy. Our study of dark matter production at the LHC goes beyond the tree-level approximation and we show examples of how higher-order corrections to dark matter production processes can affect the interpretation of the experimental results.

Abstract: Studies of dark matter lie at the interface of collider physics, astrophysics and cosmology. Constraining models featuring dark matter candidates entails the capability to provide accurate predictions for large sets of observables and compare them to a wide spectrum of data. We present a framework which, starting from a model lagrangian, allows one to consistently and systematically make predictions, as well as to confront those predictions with a multitude of experimental results. As an application, we consider a class of simplified dark matter models where a scalar mediator couples only to the top quark and a fermionic dark sector (i.e. the simplified top-philic dark matter model). We study in detail the complementarity of relic density, direct/indirect detection and collider searches in constraining the multi-dimensional model parameter space, and efficiently identify regions where individual approaches to dark matter detection provide the most stringent bounds. In the context of collider studies of dark matter, we point out the complementarity of LHC searches in probing different regions of the model parameter space with final states involving top quarks, photons, jets and/or missing energy. Our study of dark matter production at the LHC goes beyond the tree-level approximation and we show examples of how higherorder corrections to dark matter production processes can affect the interpretation of the experimental results.

Introduction
Evidence for the existence of dark matter (DM), although indirect, is quite convincing [1][2][3]. Measurements of the cosmic microwave background and baryonic acoustic oscillations predict a dominant dark matter component in the matter budget of the Universe (in the framework of standard cosmology). In addition, detection of gravitational anomalies, such as the flattening of galaxy rotation curves and the presence of gravitational lensing in the absence of visible matter (e.g. the bullet cluster [4]), strongly favours gravitational interactions of dark matter as plausible explanations.
The many hints for dark matter sparked a huge endeavour to detect it and measure its properties, leading to a number of experiments and searches which exploit very different ideas and approaches to dark matter detection. The experiments can be broadly grouped into three categories: • A wide range of underground nuclear recoil experiments aimed at detecting galactic dark matter scattering off atomic nuclei; • Searches for dark matter annihilation in the galaxy or nearby dense sources via measurements of, for instance, gamma-rays; • Collider searches in channels with large missing transverse energy ( / E T ).
However, despite an enormous experimental effort, the detection of the dark matter particles remains elusive. In fact, there is no clear indication that dark matter interacts with ordinary matter via forces other than gravity, and current experimental results are not able to put stringent bounds on the dark matter properties and couplings in a modelindependent way.
As so little is known about the true nature of dark matter, it is a useful strategy to try and constrain viable dark matter scenarios in the most model-independent way (i.e. via simplified models), confronting them with results from collider experiments, direct dark matter searches, astrophysical observations and cosmology. If or when a signal is observed, the aforementioned approach will help us to determine more accurately both the particle properties (mass, couplings, etc.) and astroparticle properties (halo properties, thermal relic density, etc.) of dark matter. Conversely, if searches result only in limits on dark matter parameters, combining constraints from different approaches aids us in excluding specific scenarios and hence narrow down the scope of viable dark matter theories.
Recent collider searches have focused mostly on studies of dark matter in the simplified model framework, where a single dark matter candidate of arbitrary spin couples to visible matter (e.g. quarks) via an s-channel or a t-channel mediator, whose quantum numbers are fixed by assumed local and global symmetries [5]. The minimal implementations of simplified dark matter models involve four basic parameters: the mass m X of the dark matter particle, the mass m Y of the mediator, the coupling constant g X of the dark matter to the mediator and the universal coupling g SM of the mediator to the visible sector (the width of the mediator is a derived quantity). Fast and efficient studies of the full simplified model parameter space require parameter scanning technology beyond simple sequential grids, due to the relatively high dimensionality of the parameter space. Past studies of simplified dark matter models have hence been limited to explorations of the parameter space in two-dimensional projections while keeping the remaining parameters fixed (see e.g. the works of refs. [6][7][8][9][10][11][12][13][14][15][16] and the references therein).
In this paper we illustrate how comprehensive studies of simplified dark matter models can be performed, exploring their full four-dimensional parameter space while taking into account constraints from collider physics, astroparticle physics and cosmology. For concreteness, we focus on a class of simplified models where the dark matter dominantly couples via a scalar mediator to top quarks (i.e. 'top-philic dark matter' scenarios). Yet, the methodology we employ is general and can be applied to other scenarios as well. We provide detailed examinations of the two-dimensional projections of the full parameter space, and we demonstrate that striking features in the structure of the viable parameter space emerge through the combination of all current constraints. We also stress that in addition to collider searches for dark matter in channels with large missing energy, in this study we also consider resonance searches in channels with fully reconstructed final states, which can be useful to constrain the properties of the mediators.
We perform the study of simplified top-philic dark matter models by using a combination of simulation tools, including the MadGraph5 aMC@NLO (MG5 aMC shorthand) event generator [17], the FeynRules package [18,19], the MadAnalysis 5 platform [20][21][22], the Delphes 3 detector simulator [23] and the MadDM program [24,25], together with an efficient parameter sampling technology based on the MultiNest algorithm [26,27]. We explore the full four-dimensional parameter space of the model in the light of existing collider and astroparticle constraints. Our analysis thus also represents a proof of concept for a unified numerical framework for comprehensive dark matter studies at the interface of collider physics, astrophysics and cosmology. This has direct implications for dark matter searches at colliders, as comprehensive phenomenological studies of dark matter models can be used to drive the experimental efforts towards the regions of the parameter space that are not already ruled out by astrophysical and cosmological constraints. In addition, we have also implemented previously unavailable experimental analyses into the MadAnalysis 5 platform, providing an added benefit of our work for future collider studies which go beyond searches for dark matter.
The article is organised as follows. Section 2 describes the details of the simplified topphilic dark matter model under consideration and discusses the constraints on the model parameter space that are implemented in our analysis setup. All cosmology and astrophysics constraints are discussed in section 3. More precisely, the relic density constraints are illustrated in section 3.1. We discuss the direct detection constraints in section 3.2, while constraints from gamma-ray flux measurements are detailed in section 3.3. Collider constraints are investigated in section 4. We study constraints from searches with and without missing transverse energy in section 4.1 and 4.2 respectively. Section 5 is then dedicated to a detailed discussion of the overall combined information coming from all the considered data. Before concluding in section 7, we briefly elaborate in section 6 on whether the potential diphoton excess observed by the ATLAS and CMS collaborations at a diphoton invariant mass of m γγ 750 GeV [28,29] could be interpreted within the considered class of simplified top-philic dark matter models. We provide more information on the mediator width in appendix A. As a validation of our calculations, we perform a detailed comparison between MadDM and MicrOMEGAs in appendix B.1, give details on the annihilation cross section of dark matter in the top-philic model in appendix B.2 and present the validation of the CMS tt + / E T and monojet implementation in the Mad-Analysis 5 framework in appendix C.

Simplified top-philic dark matter model and its numerical implementation
The simplified top-philic dark matter model that we consider is constructed by supplementing the Standard Model (SM) with a Dirac-type fermionic dark matter candidate X and a scalar mediator Y 0 . The interactions of the two particles are described by the Lagrangian where the new physics interaction strengths are denoted by g t and g X for the mediator couplings to the Standard Model sector and to dark matter respectively. We have assumed an ultraviolet-complete description of the scalar theory where the mediator couples to quarks with a strength proportional to the Standard Model Yukawa couplings, so that we neglect all light quark flavour couplings and only include the coupling of the mediator to the top quark, y t = √ 2m t /v where v = 246 GeV is the Higgs vacuum expectation value and m t is the top quark mass. Note that the model in eq. (2.1) is neither complete, nor stable under radiative corrections. Couplings to the top quark induce a mixing with the standard model Higgs, which we set to zero by construction. In addition, loop corrections will also generate finite couplings to pairs of W and Z bosons, which we will omit in the following.
The model contains four free parameters (two couplings and two masses), while the width Γ Y is fixed by the remaining model parameters. In addition to the Lagrangian of eq. (2.1), we could also have considered mediator couplings to leptons. They however cannot be well constrained by LHC searches and dark matter direct detection data, and we have excluded them from our model description. We will nonetheless comment on their relevance for relic density predictions and dark matter indirect detection signals in sections 3.1 and 3.3. The Lagrangian of eq. (2.1) induces dimension-five couplings of the mediator to gluons and photons via loop diagrams of top quarks. The loop-induced operators can be relevant in the context of both astrophysical and collider searches for dark matter. The couplings of the mediator to gluons and photons are given, at the leading order (LO), by the effective operators with the effective couplings being In the above expressions, Q 2 denotes the virtuality of the s-channel resonance, while F S is the one-loop form factor with F S (x) → 2/3 for x 1. Eq. (2.4) contrasts with the Standard Model Higgs case where the effective Higgs-photon coupling receives contributions from vector-boson loopdiagrams that are absent in our simplified dark matter model setup. As a result, the gluon and photon effective couplings to Y 0 are characterised by a larger hierarchy compared to their Higgs counterparts.
The tree-level partial decay widths of the scalar mediator are given by where β t,X = 1 − 4m 2 t,X /m 2 Y and Θ(x) is the Heaviside step function, and we ignored the top quark width in the expression for Γ(Y 0 → tt). The loop-induced Y 0 partial widths are The Y 0 partial width to photons is by construction always smaller than the partial decay width into a pair of gluons by virtue of α 2 s /α 2 e ∼ 100. In addition to a coupling suppression, other decay processes such as the loop-induced Y 0 decays into Zγ, ZZ and hh final states receive a kinematic suppression. Couplings of Y 0 to ZZ and hh could also appear at tree level in our model, but in the spirit of simplified models, we define them to be vanishing. In the following we hence safely approximate the total decay width for the mediator to be the sum of eqs. (2.6), (2.7) and (2.8).
The total decay width and the branching ratios of the mediator into tt, XX, gg and γγ final states are displayed in figure 1 for different choices of new physics couplings and masses. Light mediators with masses below the top-quark pair or the dark matter pair decay thresholds are narrow states, while above these thresholds, large Γ Y /m Y values are possible in particular for large couplings. For mediators with m Y m t , m X , the dominant decay channel is into a pair of gluons. In contrast, heavy mediators with mass m Y > m t , m X decay predominantly into pairs of top quarks and/or dark matter particles, where the exact details of the partial width values strongly depend on the masses and couplings. The branching ratio of Y 0 to photons is always suppressed, as argued above. We present in appendix A the dependence of the Γ Y /m Y ratio on the g t and g X couplings for different mass choices and on the m Y and m X masses for different coupling choices.
Our top-philic dark matter model can be probed in different ways including astrophysical and collider searches, as listed in table 1. The relative importance of the various searches depends on the hierarchy of the dark matter, mediator and top-quark masses, as well as on the hierarchy between the couplings. Starting with the dark matter relic density, the annihilation cross section is dominated by subprocesses with top-quark final states for m X > m t , and by annihilation into gluons and to a lesser extent photons for light dark  Figure 1. Ratio of the mediator width to its mass Γ Y /m Y (upper panels) and mediator branching ratios (lower panels) as a function of the mediator mass for different coupling choices and a dark matter mass fixed to m X = 50 GeV (solid lines) and 300 GeV (dashed lines). matter particles with m X < m t . If the mediator is lighter than the dark matter state, an additional annihilation channel into a pair of mediators can open up. The annihilation mechanisms into top-quarks, gluons/photons and mediators moreover provide an opportunity to indirectly search for dark matter, e.g. in gamma-ray data. The interactions of the dark matter particles with nuclei, relevant for direct detection experiments, proceed via mediator exchanges. The mediator-nucleon coupling is in turn dominated by the scattering off gluons through top-quark loops. Dark matter production at the LHC proceeds either through the production of the mediator in association with top quarks, or from gluon-fusion through top-quark loops. Searches at the LHC can be classified into two categories regarding the nature of the final states that can contain missing transverse energy / E T or not. Searches involving missing energy may include final state systems containing a top-quark pair and probe in this way the associated production of a top-antitop-mediator system where the mediator subsequently decays into a pair of dark matter particles. Alternatively, the mediator can be produced via gluon fusion through top-quark loops, where the probe of the associated events consists of tagging an extra radiated object. This yields the well-known monojet, mono-Z and mono-Higgs signatures. We do not consider the monophoton channel, as photon emission is forbidden at LO in our simplified model by means of charge conjugation invariance. The second search category is related to final states without any missing energy, i.e. when the mediator decays back into Standard Model particles. This includes decays into top-quarks, leading to final states comprised of four top quarks, into a top-quark pair, as well as into a dijet or a diphoton system via a loop-induced decay. This is, however, relevant only for on-shell (or close to on-shell) mediator production.
We proceed with a description of the numerical setup for our calculations. In the following sections, we explore the full four-dimensional model parameter space and present results in terms of two-dimensional projections. We perform the four-dimensional sampling using the MultiNest algorithm [26,27], where we assume Jeffeys' prior on all the free parameters in order not to favour a particular mass or coupling scale. The choice of prior ranges for the parameters is summarised in table 2, in which we have chosen to limit the coupling values to a maximum of π to ensure perturbativity. We implement the relic density constraints into MultiNest using a Gaussian likelihood profile, while for the direct detection limits we assume a step likelihood function smoothed with half a Gaussian. In addition, the sampling imposes that the model is consistent with values of Γ Y such that the mediator Y 0 decays promptly within the LHC detectors. Table 3 summarises the constraints that we have imposed on the model parameter space.
Throughout our study, we assume that X is the dominant dark matter component, namely that it fully accommodates a relic density Ω DM h 2 as measured by the Planck satellite [30]. Concerning the direct detection of dark matter, we consider the currently most stringent bounds on the spin-independent (SI) nucleon-DM cross section as measured by LUX for dark matter with m X > 8 GeV [31] and by CDMSLite for 1 GeV< m X < 8 GeV [32]. In section 3.3, we focus on indirect detection constraints that are imposed on the basis of the gamma-ray measurements achieved by the Fermi-LAT telescope [33,34]. Those bounds are however not applied at the level of the likelihood function encoded in our MultiNest MultiNest parameter  We derive collider constraints on the simplified top-philic dark matter model using the MG5 aMC [17] framework and the recast functionalities of MadAnalysis 5 [20][21][22] (where appropriate). We apply the LHC constraints on the top-philic dark matter model with two different procedures. On one side, similarly to what has been performed for the indirect detection bounds, we reprocess the scenarios that accomodate the observed relic density and that are compatible with LUX and CDMSLite data. However, we also study the collider bounds on the parameter space independently of any astrophysics and cosmology consideration and by relaxing the narrow width requirement (allowing Γ Y /m Y to be of O(1)) as well. In order to increase the sensitivity of the LHC searches, we allow for wider coupling ranges of 10 −2 < g X < 2π and 10 −2 < g t < 2π. The collider study without any cosmological and astrophysical constraint therefore includes the cases where the dark matter is not a standard thermal relic (i.e. its relic density is a result of a nonthermal mechanism or a non-standard evolution of the Universe). Details are provided in section 4 and appendix C for what concerns the validation of the CMS analyses that we have implemented in MadAnalysis 5 for this work.
In conclusion to this section, we point out that even though our current work focuses on a dark matter candidate which is a Dirac fermion, a more general implementation of simplified dark matter models in FeynRules [18,19] can also account for pseudoscalar mediators as well as for CP -mixed states and for dark matter particles which are real or complex scalars [36][37][38]. The corresponding model files have been used in this work and can be downloaded from the FeynRules model repository [39] that also includes a model where the mediator is a spin-1 state that couples to either a fermionic or a scalar dark matter candidate [36]. All the models allow for the automated calculation of next-toleading-order (NLO) effects and loop-induced leading-order (LO) processes in QCD in the context of LHC predictions.

Cosmological and astrophysical constraints
We begin our analysis of the simplified top-philic dark matter model with a detailed discussion of the cosmological and astrophysical constraints.

Constraints from dark matter relic density
Dark matter annihilation in the early Universe is determined, in the simplified top-philic dark matter model, by a combination of three processes, where we have omitted the annihilation into photons as it is always suppressed compared to the annihilation into gluons. The analytic expressions for the thermally averaged annihilation cross section in the non-resonant region σv rel corresponding to each of the processes listed above are provided in appendix B.2. The first two processes proceed via an s-channel Y 0 exchange (first two rows of table 1), while the third process consists of a t-channel X exchange (third row of table 1). The resonance structure of the s-channel processes implies that the width of Y 0 potentially plays an important role in the determination of the relic density assuming a dominant annihilation via the processes (I) and (II), while the effects of the Y 0 width are mostly negligible if the annihilation dominantly proceeds via the t-channel X exchange process (III). According to the hierarchy between the dark matter mass m X , the mediator mass m Y and the top quark mass m t , different situations can occur. Qualitatively, one expects that: • for m Y m X m t : process (I) is dominant as the tree-level annihilation into a pair of top quarks is kinematically allowed, the annihilation into gluons being loop suppressed, and the one into a pair of mediators kinematically suppressed; • for m X m t , m Y : dark matter annihilates into a pair of gluons as in process (II), since it is the only kinematically allowed channel; • for m t m X m Y : relic density is determined by process (III) since annihilation into top quarks is kinematically forbidden and the one into gluons occurs away from the resonant pole of m Y ; • for m X > m t , m Y and m Y < 2m t : similarly to the case above, the dominant annihilation mechanism is process (III), as annihilation into top quarks occurs far from the resonant pole and is suppressed kinematically; • for m X > m t , m Y and m Y > 2m t : processes (I) and (III) are competitive and the dominant process among the two is determined by the hierarchy between the g t and g X couplings.
Requiring our simplified top-philic dark matter model to result in a dark matter relic density consistent with the most recent Planck measurements [30] implies strong constraints on the viable regions of the parameter space. As an illustration, we consider the region of the parameter space in which m t m X m Y , where we expect the dominant annihilation mechanism of dark matter to be process (III) and to give rise to a pair of mediators. In this region, the thermally averaged annihilation cross section approximately reads so that it is clear that imposing that the relic density predictions agree with Planck data leads to a stringent constraint on the ratio g 2 X /m X . The argument is more involved in parameter space regions where the total mediator width Γ Y plays a role, as the relevant quantity involved in the relic density calculation is in general not σv rel ann but dx σv rel ann (x) where x ≡ m X /T and σv rel ann is a non trivial function of x. This is especially true, for instance, for the Breit-Wigner-type amplitudes that appear in processes (I) and (II).
In order to provide a more detailed quantitative analysis, we have performed a fourdimensional scan the top-philic dark matter model parameter space and examined the effects of imposing relic density constraints on the allowed/ruled out parameter sets. Figure 2 reveals the rich structure of the four-dimensional parameter space allowed by relic density measurements. The bulk of the allowed parameter points lies in the region where m X > m Y , and the annihilation cross section is dominantly driven by process (III). This region of the parameter space has the particularity of not being reachable by traditional monojet, monophoton, mono-Z and mono-Higgs searches at colliders. The decay of the mediator into a pair of dark matter particles is indeed not kinematically allowed, so that any new physics signal will not contain a large amount of missing energy. The model can however be probed at colliders via dijet, diphoton, tt (plus jets) and four-top analyses. We elaborate on this point more in section 4.2. The characteristic mediator width Γ Y in this region tends to be extremely small, with values of at most 10 −4 GeV as shown in the top left panels of figure 2. This is expected as the width is mostly controlled by the decays into gluons, and into top quarks in the regions where this decay is kinematically allowed, the decay into a pair of dark matter particles being forbidden.
In the region where m X m t and m Y 2m t , the mediator decay into a tt final state is kinematically allowed and the dark matter annihilation cross section is driven by the XX → Y 0 → tt process. The only other parameter space region that is not ruled out by the relic density data is centered around the resonance region where m Y ∼ 2m X . The extension of the region away from the resonance pole is due to the Y 0 width that can reach O(10) GeV. The resonant region extends to lower m X and m Y values, and is the only allowed region when both m X and m Y are smaller than m t (but with Y 0 decays into a pair of dark matter particles being allowed). This has interesting implications for LHC searches as the low dark matter/mediator mass region is the one where colliders have the best sensitivity, in particular through monojet searches (see for instance section 4.1.2). Relic density constraints favour g X couplings of O(1) in most of the scanned parameter range as evident in the second and third panels of figure 2, regardless of the actual value of the g t coupling which is irrelevant in the m X m Y region (upper right panel of figure 2) as it does not enter the calculation of the relic density.
The structure of the ruled out parameter space regions shows several other interesting qualities. The most striking feature is that almost the entire region where m Y 2m X does not lead to predictions of a dark matter relic density in agreement with the observations. There are also no allowed points for m X m t , except very close to the resonance line. . Results of our four-dimensional parameter scan using MadDM, projected onto the (g X , g t ) plane (left) and (m X , g t ) plane (right). All represented points feature a relic density in agreement with Planck data and a narrow width mediator.
This region is characterised by a dominant mediator decay into gluons, which results in typical Γ Y /m Y 1, a small total dark matter annihilation cross section, and hence an overproduced dark matter. The upper limit imposed on the size of the couplings (see table 2) is largely responsible for the absence of allowed points in the region. For instance, taking any m X value so that the predicted relic density agrees with the observed value, an increase in m Y will result in a decrease of the annihilation cross section, in turn leading to a higher relic density. The only way (away from the resonance) to restore the correct relic density is then to increase the size of g X and/or g t . However, our results show that even for couplings of O(1), the cross section in this region is too small not to overproduce dark matter.
The region of parameter space between m Y ∼ m X and m Y ∼ 2m X is consistent with the above-mentioned argument. This strip of the ruled out parameter space can be seen as a part of the larger ruled out region for which m Y m X . Tuning m X to be close to m Y /2 and assuming a relatively small Γ Y value is the only way to enhance the dark matter annihilation cross section and not overproduce dark matter.
In addition to projections of the allowed parameter space onto the (m Y , m X ) plane, we have also studied several other projections. Figure 3 shows the projections of our results onto the (g X , g t ) plane (left) and (m X , g t ) plane (right), where we show m X and m Y as a colourmap in the first and second panel respectively. Regardless of the value of m X and m Y in the considered scan range, there are no solutions for g X and g t which satisfy the relic density constraint in the region where g t 10 −2 and g X 10 −1 . This finding is consistent with the left lower panels of figure 2 where we have found that a correct relic density favours g X couplings of O(1). Furthermore, we can observe that only couplings of GeV, while in the majority of the allowed (g X , g t ) parameter space regions Γ Y /m Y 1. We find no striking features in the (m X , g t ) projection of the scanned parameter space.
The unpopulated regions in the lower left corners are artifacts of the lower limit on the coupling size of 10 −4 .
As a validation, we have cross checked our calculations with the MicrOMEGAs code. The results obtained with MadDM and MicrOMEGAs agree in most of the parameter space, except in the region where g t and m X are small. Some numerical discrepancies are expected to occur in this region, as shown in appendix B.1 and by comparing figures 2 and 20.
As a last remark, allowing the scalar mediator to couple to all quarks and leptons would only have a minor impact on our results. The region dominated by the process (III) will indeed stay unchanged, since it is insensitive to the coupling between Y 0 and the Standard Model fermions. As far as it concerns dark matter annihilation via an s-channel Y 0 exchange, one would have to sum up over all the possible final states kinematically open. This would increase the total annihilation cross-section and decrease Ω DM h 2 , implying that the constraint of having Ω DM h 2 ∼ 0.12 leads to a rescaling of all fermionic couplings towards smaller values with respect to the g t values shown in this work. The major difference would reside in a potentially larger decay width for Y 0 and hence wider "bands" of the resonance regions of the allowed parameter space.

Constraints from direct detection
Simplified models of dark matter which feature couplings to quarks and gluons can also be bounded by results from underground direct dark matter detection experiments. In topphilic dark matter scenarios, dark matter scatters off nucleons via the t-channel exchange of Y 0 , where the scattering off gluons via triangle top loops accounts for the dominant contribution to the DM-nucleon scattering rate.
The spin independent (SI) dark matter-nucleon cross section is given by [40] 1 is the gluon form factor and the sum runs over the light quarks q = u, d, s, where m n ≈ 0.938 GeV is the nucleon mass and m t = 173 GeV is the top quark mass. The expression in eq. (3.2) does not depend on Γ Y , simplifying the constraints which can possibly be derived from direct detection. For instance, considering a scenario in which generic m X and m Y masses are fixed and where the dominant annihilation process is process (I), direct detection directly constraints the product g X g t . Extracting the constraint on this quantity in a generic fashion is much more complicated in the case of dark matter annihilation in the early Universe and at colliders, as the processes involved in dark matter relic density and dark matter production calculations intrinsically depend on a quantity which is proportional to g X g t /Γ Y . The running of the g X and g t couplings could have an effect on the value of the spin independent DM-nucleon scattering cross section [41,42]. However, a proper inclusion  of the running couplings would require a careful treatment of the renormalisation group evolution via multiple energy scales which is beyond the scope of our current effort. Instead, we restrict here our calculations to constant g X and g t values. The effect of the running couplings would then be equivalent to a rescaling of g X and g t to different values. Next, we have repeated the four-dimensional parameter scan from section 3.1 including into the MultiNest likelihood function also bounds stemming from direct detection. Figure 4 shows the results of the scan, projected onto three different planes, where we removed the points excluded by the 95% confidence limit (CL) bound from LUX and CDMSLite. Direct detection rules out a major portion of (m Y , m X ) space allowed by the relic density constraints (regardless of the coupling value) in the region where m X m Y , where collider bounds are irrelevant. Figure 4 hence serves as a good example for the complementarity among direct detection, relic density and collider bounds. In the (g X , g t ) plane, direct detection does not rule out a well defined-region (top-right panel of figure 4), indicating that for any pair of couplings (g X , g t ) in the range of [10 −4 , π] which are allowed by the relic density constraint, it is always possible to find a pair of (m Y , m X ) values which are not ruled out by direct detection data. In the (m X , g t ) projection, we finally observe that direct detection rules out a well-defined portion of the parameter space. Furthermore, the constraint also rules out small width points for g X 0.1 and m X m t . Direct detection bounds are indeed more sensitive to dark matter masses in the ballpark of 10 to 200 GeV and quickly deteriorate at larger dark matter masses, since the event rate in the detector scales as 1/m 2 X . We also see that the direct detection exclusion limit is able to rule out a large portion of the parameter space where Y 0 is light, below 30 GeV, while the sensitivity is quickly lost for heavier masses of the scalar mediator. This can be understood by the 1/m 2 Y dependence of the SI elastic cross section of eq. (3.2). Both mass dependences are illustrated by the lower panel of figure 4.

Constraints from indirect detection
Top-philic dark matter annihilation in the present Universe could result in fluxes of cosmic rays and prompt gamma-rays, which can also be used to infer useful limits on the model parameter space. The annihilation of a XX pair in the galactic halo (or in dense environments of galactic centers) and the subsequent production of a secondary gamma ray flux is dictated by the same processes (I), (II) and (III) that set the relic abundance. These processes give rise to a continuum of secondary photons due to the decay and subsequent QED showering of the pair-produced top quarks, gluons and/or mediators. As already mentioned in section 2, a direct coupling of the mediator to a pair of prompt photons is induced at higher order in perturbation theory via a loop of top quarks. Hence, analogously to process (II), the process XX → γγ exists and yields the production of two monochromatic photons that could be detected in searches for lines in the gamma-ray spectrum. 2 Finally, photons arising from process (III) and the subsequent decay of the mediator into two photons do not provide a signal line as the mediators are in general not produced at rest in the annihilation process.
Similarly to the relic density case, measurements of the gamma-ray fluxes can potentially constrain the coupling g X for the t-channel process (III) or the product of couplings g X g t in the case of an s-channel annihilation via the processes (I) and (II). However, it is important to highlight the differences between factors which are constrained by the dark matter relic density and by its indirect detection. The relic density is an integrated result over the thermal history of the Universe. Hence, the width of the resonance is important, Conversely, the characteristic velocity of the dark matter particles today is of the order of v ∼ 10 −3 , implying highly non-relativistic dark matter annihilation. The width of the mediator in an s-channel dark matter annihilation process is hence relevant for indirect detection only in the case of . Dark matter annihilation cross section at present time that is relevant for gamma-ray limits extracted from dwarf spheroidal galaxies measurements (left) and gamma-ray line searches (right). We show a maximal estimate of (σv rel ) tot and (σv rel ) γγ obtained by where v ∞ is the escape velocity for dwarf spheroidal galaxies and the galactic center, respectively. All represented points are compatible with the relic density, a narrow width mediator and the direct detection requirements.
Searches for gamma-ray signals of dark matter annihilation weakly constrain our simplified top-philic dark matter model. We have investigated results from gamma-ray line searches in the inner galactic region [34], as well as continuum gamma-ray measurements from dwarf spheroidal galaxies [33] and found no meaningful exclusion of the parameter space once the relic density and direct detection constraints are imposed. The lack of additional useful bounds is expected, as the annihilation of dark matter in the present Universe is p-wave suppressed, i.e. σv rel ∝ v 2 rel for all three annihilation channels (see appendix B.2 for more detail). This contrasts with scenarios in which the mediator is a pseudoscalar state that implies that the p-wave suppression at low dark matter velocity is only present for process (III), so that the gamma-ray constraints should be significantly stronger.
The gamma-ray line searches constrain the velocity-averaged cross section for the direct dark matter annihilation into two photons. Due to its p-wave suppression, this quantity is very sensitive to the choice of the velocity distribution of the dark matter in the galaxy which is subject to large uncertainties (see e.g. ref. [43]). We adopt a conservative viewpoint here, evaluating the annihilation cross section at the highest possible velocity v rel = 2v ∞ with v ∞ being the escape dark matter velocity for our galaxy which we take to be v ∞ = 550 km/s [44]. The left panel of figure 5 shows the respective result for (σv rel ) γγ . The limits from gamma-ray line searches lie between 2 × 10 −32 cm 3 s −1 (for dark matter masses around 1 GeV) and 4 × 10 −28 cm 3 s −1 (for dark matter masses around 500 GeV).
Searches for gamma-ray signals in dwarf spheroidal galaxies constrain the total the annihilation cross section at (two times) the escape velocity, the escape velocity of the considered dwarf spheroidal galaxies being typically much smaller and of the order of 10 km/s [45], which leads to a heavy suppression of the dark matter annihilation cross section. The right panel of figure 5 shows the annihilation cross section evaluated for v ∞ = 50 km/s. The cross sections are much smaller than the constraints which are around 10 −26 cm 3 s −1 , the exact details depending on the dark matter mass and the relevant annihilation processes.
In cases where we would have allowed for leptonic couplings of the scalar mediator Y 0 , our general conclusion about the poor ability of indirect dark matter searches to constrain the model parameter space remains unchanged. Dark matter annihilation into leptonic final states could give rise to additional continuum gamma-ray or positron fluxes, but the overall normalisation of σv rel would not change significantly and remain four to five orders of magnitude below the current bounds. Even under the most aggressive assumptions, all obtained bounds would still be far from being able to constrain a top-philic dark matter model with scalar mediators.

Collider constraints
As discussed in section 2, simplified top-philic dark matter scenarios can be probed at colliders through the production of the mediator either in association with a top-quark pair or through a top-quark loop. Depending on the mass and coupling hierarchy, the mediator decays either into a pair of dark matter particles, which results in signatures including missing transverse energy ( / E T ), or into Standard Model final states. The size of the cross sections associated with these two classes of mediator production mechanisms is depicted in figure 6 where we present their dependence on the mediator and dark matter masses m Y and m X . For the case where the mediator is singly produced, we use the Higgs cross section values that are reported in the Higgs Cross Section Working Group documentation [46] and that are evaluated at the next-to-next-to-leading order (NNLO) accuracy in QCD. For all the other cases, the hard-scattering cross section is convoluted with the NNPDF 2.3 [47] set of parton distribution functions (PDF) within MG5 aMC, the PDFs being accessed via the LHAPDF library [48,49]. We employ a five-flavournumber scheme, and leading-order (LO) and next-to-leading-order (NLO) PDFs are used where relevant. The renormalisation and factorisation scales are set to half the sum of the transverse mass of all the final-state particles both for LO and NLO calculations, and the scale uncertainty is estimated by varying the two scales independently by a factor of two up and down. Additional details on the calculation of the Y 0 tt cross section are provided in ref. [36] while loop-induced processes are extensively documented in ref. [37].
All the cross sections shown in figure 6 are proportional to g 2 t and we therefore arbitrarily choose g t = 1 as a benchmark. In this case, sizeable cross sections of 10 1 − 10 3 pb are expected for the production of light mediators with m Y 100 GeV at a centre-of-mass energy of 8 TeV (left panel), the dominant mechanism being the loop-induced gg → Y 0 production mode. Requiring an extra hard jet in the final state reduces the cross section by a factor which depends on the missing energy (or the jet transverse momentum p T ) selection, and the production rates are not sensitive to the mediator mass as soon as the latter is smaller than the / E T selection threshold. The cross sections for producing the mediator in association with a Standard Model Higgs or Z boson are further suppressed. In contrast, the cross section related to the production of the mediator in association with a top-quark pair is significant for light mediators, but falls off quickly with the increase in  the mediator mass due to phase-space suppression. As a result, a change in the collider energy from 8 to 13 TeV is important for heavy mediators and the cross section can be enhanced by about an order of magnitude. In the right panel of figure 6, we further show first that the cross sections are constant when the dark matter particle pair is produced through the decay of an on-shell mediator, and next that they are considerably suppressed when the mediator is off-shell, especially for the ttXX channel. As already mentioned, the collider searches which provide the most relevant constraints on simplified top-philic dark matter models are based on the production channels shown in figure 6 and can in general be divided into two categories. The first category involves signals with missing transverse energy originating from the production of dark matter particles that do not leave any trace in the detectors and that are accompanied by one or more Standard Model states. The most relevant searches of this type are the production of dark matter in association with a top-quark pair and the loop-induced production of dark matter in association with a jet, a Z boson or a Higgs boson. This is discussed in section 4.1. The second category of searches relies on Y 0 resonant contributions to Standard Model processes. In our scenario, dijet, diphoton, top-pair and four-top searches are expected to set constraints on the model parameter space. This is discussed in section 4.2. As shown below, missing-energy-based searches and resonance searches are complementary and necessary for the best exploration of the model parameter space at colliders.  [53] h → bb decay jj σ(m Y = 500 GeV) < 10 pb CMS [54] Only when m Y > 500 GeV γγ σ(m Y = 150 GeV) < 30 fb CMS [55] Only when m Y > 150 GeV tt σ(m Y = 400 GeV) < 3 pb ATLAS [56] Only when m Y > 400 GeV tttt σ < 32 fb CMS [57] Upper limit on the SM cross section Table 4. Summary of the 8 TeV LHC constraints used in this paper.
In the rest of this section, we study collider constraints independently from the cosmological and astrophysical ones, and we dedicate section 5 to their combination. We moreover allow the mediator couplings to be as large as 2π and do not impose any constraint on the mediator width over mass ratio. We summarise the relevant 8 TeV LHC constraints used in this study in table 4 and give details on the tt + / E T and monojet searches that have been recast in the MadAnalysis 5 framework in appendix C. Dark matter production in association with a top-quark pair (tt + / E T ) has been explored by both the ATLAS [58] and CMS [59] collaborations within the 8 TeV LHC dataset, and limits have been derived in particular in the effective field theory approach [60,61]. Such analyses could however be used to derive constraints in other theoretical contexts, and we choose to recast the CMS search to constrain the parameters of the simplified top-philic dark matter model under consideration. In this work, we simulate ttXX events at the NLO accuracy in QCD by making use of MG5 aMC. The first study of the genuine NLO effects on the production of a system composed of a pair of top quarks and a pair of dark matter particles has been presented in ref. [36] in which NLO K-factors have been investigated both at the total cross-section and differential distribution level for a series of representative benchmark scenarios. Here, we explore the impact of the NLO corrections on the exclusion limits originating from the tt + / E T channel. In order to examine the reach of the CMS search, we start by performing a twodimensional scan of the mediator and dark matter masses with fixed mediator couplings, similar to figure 7 in ref. [7]. The same scan is performed at both LO and NLO accuracy concerning the simulation of the hard scattering process, which allows us to determine the impact of the QCD corrections on the exclusion bounds. Before presenting the results for the excluded regions and to facilitate the discussion, we show the dependence of the LO cross section for g t = g X = 4 on the new physics masses and the corresponding K-factors in figure 7. The cross section is the largest in the low mass regions where the mediator can resonantly decay to a pair of dark matter particles, and falls steeply in the off-shell regions. In particular, the region where 2m X < m Y < 2m t is characterised by mediator decays either into a pair of dark matter particles or into a pair of gluons. These two decay rates are related by (see section 2) which suggests that the decay rate into a pair of dark matter particles is always significantly higher, except in the case of a large hierarchy between the couplings (g t /g X 100). For m Y > 2m t , the Y 0 → tt decay mode is open, and tt + / E T production turns out to be suppressed by the visible decay channels of the mediator, unless g X > g t . Such a feature has already been illustrated in figure 1. The NLO K-factors related to tt + / E T production (right panel of figure 7) are found to vary from 0.96 to 1.15 in the range of masses examined here, the QCD corrections being more important in the low mass region.
The results for the exclusion regions are shown in figure 8 when LO (left panel) and NLO (right panel) simulations are used; see more details on the recasting procedure in appendix C.1. Setups excluded at the 40%, 68% and 95% confidence level (CL) are marked separately in the figures. As expected from the total cross section results, all excluded points (at the 95% CL) lie in the triangular low-mass region where the mediator resonantly decays into a dark matter particle pair. The exclusion region reaches mediator masses of about 200-250 GeV if close to threshold (m Y ∼ 2m X ). This region is in fact not exactly triangular as for a given mediator mass, not all dark matter masses below m Y /2 are excluded. This is related to the parametric choice of g t = g X = 4 for which the mediator width can become large. In this case, the narrow width approximation is not valid and the tt + / E T cross section acquires a dependence on the dark matter mass even in the resonant region.
Comparing the LO and NLO results, we observe that in the low mass resonant region where the K-factor is small and of about 1.10, the exclusion contours are mildly modified and this small 10% shift in the cross section does not lead to any significant change. For larger mediator masses, the K-factors are ∼ 1 and therefore do not imply a modification of the exclusion regions, if the central prediction at the default choice of scale is considered. However, the inclusion of NLO corrections significantly reduces the theoretical error and thus leads to sharper exclusion bounds as discussed below.
In order to further investigate the effects of the NLO corrections, we select three benchmark scenarios for which we perform a detailed study. These benchmarks are defined in table 5 where they are presented along with the corresponding LO and NLO cross-sections and the CL exclusion obtained with MadAnalysis 5. As discussed in appendix C.1, the most relevant observables for this analysis consist of the / E T , M T ( , / E T ) and M W T 2 for which distributions are shown in figure 9. We normalise the distributions to 100, 10 and 1 for the scenarios I, II and III respectively to ensure that they are all clearly visible in the figure. Moreover, we also indicate the scale uncertainty bands that have been obtained from a scale variation of 0.5µ 0 < µ R,F < 2µ 0 . In agreement with the findings of ref. [36], higher-order corrections have a rather mild effect on the distribution shapes for all key observables. Using NLO predictions however leads to a significant reduction of the scale uncertainties compared to the LO case. In table 5, one can also see that the use of NLO predictions leads to a significant reduction of the uncertainty in the cross section which propagates down to the CLs. NLO predictions therefore allow us to draw more reliable conclusions on whether a parameter point is excluded.

Mono-X final states
In addition to the constraints that can be derived by means of tt+ / E T probes and that have been discussed in the previous section, mono-X searches can also be relevant for obtaining bounds on our top-philic dark matter model. Monojet [51,62,63] Table 5. Benchmark scenarios used to investigate the impact of the NLO corrections on the tt + / E T CMS search. The LO and NLO cross sections at 8 TeV LHC are shown together with the CL exclusion obtained from MadAnalysis 5. The uncertainties originating from scale variation (0.5µ 0 < µ R,F < 2µ 0 ) are also shown. this work. In contrast to tree-level dark matter production in association with a pair of top quarks, the production of a pair of dark matter particles with a jet, a Z-boson or a Higgs boson proceeds via a gluon fusion top-quark loop diagram. Although they have been largely studied by ATLAS and CMS, monophoton analyses cannot be used as charge conjugation invariance forbids the existence of a monophoton signal for the spin-0 mediator scenario.

Monojet
We start by discussing constraints that can be imposed by the CMS 8 TeV monojet analysis [51]. For this study, hard-scattering events are generated at the LO accuracy within MG5 aMC, and the matching with parton showers is made with Pythia 6. The results are analysed in MadAnalysis 5 that also takes care of the detector simulation using its interface with Delphes 3. This recasting procedure allows us to exclude any specific parameter space point at any desired confidence level, our exclusion being conservatively derived on the basis of the signal region that drives the strongest bound. This limitation is related to the lack of public information, the statistical model used by CMS for the combination being not available. One can find more details for the recasting procedure in appendix C.2.
Similar to the tt + / E T analysis of the previous section, we perform a two-dimensional scan on the mediator and dark matter masses while fixing both new physics couplings to g t = g X = 4 (as in figure 5 in ref. [7]). Figure 10 shows our results, where we represent the scenarios excluded at the 40%, 68% and 95% CL. The bulk of the excluded points lie again in the triangular low-mass region where the mediator resonantly decays into a pair of dark matter particles. Except for the small subset of points excluded at the 40% and 68% CL in the region where m Y < 2m X , the extent of the exclusion region is determined by the significant reduction of the monojet cross section below the resonant production threshold already presented in figure 6. The pp → Y 0 j cross section indeed rapidly falls with m Y , reaching levels beyond the sensitivity of the 8 TeV search at m Y ∼ 500 GeV. In addition to the decrease of the Y 0 j production cross section, the opening of the mediator decay mode into a top-antitop system when m Y > 2m t leads to a further reduction of the monojet production rate. In comparison with the tt + / E T case, the monojet search overall    Table 6. Benchmarks used to investigate the differential distributions related to the CMS monojet analysis. The corresponding cross sections for a / E T > 150 GeV selection are shown in the second column.
appears to be more constraining, especially for higher mediator mass values thanks to the larger monojet cross section.
As shown in ref. [37], the shape of key monojet differential distributions differs in the resonant and in the off-shell parameter space regions. While the total cross section falls dramatically in the off-shell region m Y < 2m X (as shown in figure 6), the / E T and jet transverse momentum distributions tend to be harder for off-shell production. We demonstrate this feature with a detailed investigation of three benchmark points defined in table 6. They consist of two resonant scenarios with different mediator masses and one non-resonant scenario. The monojet production rate is also indicated in the table, and we present normalised distributions relevant for the monojet analysis in figure 11. The off-shell scenario yields harder distributions compared to the resonant cases. This implies that a larger fraction of events features high missing transverse energy ( / E T >250 GeV) and populates the different signal regions of the CMS analysis. As a result, a better sensitivity is found than what one might expect from considering the total cross section alone. This feature leads to the exclusion of dark matter scenarios where m Y < 2m X , as depicted in figure 10.
In our simulation of the monojet signal, we have ignored the possible impact of the merging of event samples featuring different final state jet multiplicities. A reliable de- scription of the high transverse momentum spectra of the leading jet typically necessitates the merging of event samples including at least one and two jets in the final state [37]. We have explicitly verified that for both resonant and off-shell scenarios, employing a merged sample does not have a big impact on the / E T distribution and therefore on the resulting exclusion contours. This originates from the analysis selection strategy that requires one single hard jet and rather loose requirements on the second jet, so that the configuration that dominates consists of a single hard jet recoiling against the missing energy. Such a configuration is described similarly by the one-jet and merged samples. We nevertheless stress that the importance of the merging procedure has to be checked on a case-by-case basis as this depends on the analysis, so that higher multiplicity samples might be necessary to accurately describe the relevant distributions.

Mono-Z and mono-Higgs
In addition to the use of monojet processes, we explore the possibility of constraining the parameter space of our model using mono-Z and mono-Higgs production. While the production rates are much smaller than the monojet rate as seen in figure 6, the backgrounds can be also small. Therefore, these search channels can be sensitive to the top-philic simplified dark matter model, as we will see below. Here, instead of employing a full recasting procedure as in the tt + / E T and monojet analyses, we perform parton-level analyses to provide rough estimates of the constraints on our model parameters.
We rely on the CMS search for dark matter production in association with a Z-boson that decays leptonically [52], in which a 95% CL upper limit on the visible cross section of 0.85 fb is obtained once a / E T requirement of at least 150 GeV and the minimal detector selection requirements for the leptons (p T > 20 GeV and |η | < 2.5) are considered. We generate events for this process, and after applying the above fiducial selection requirements we obtain a cross section of 0.30 fb for (m Y , m X ) = (100, 10) GeV and g t = g X = 1. We show in figure 12 the / E T and leading lepton transverse momentum distributions (red lines) without and with applying the selection strategy. While we have not performed a detailed study, simple estimates show good prospects for setting limits on the parameter space of the model using the mono-Z analysis results. Using the upper limit of 0.85 fb, scenarios with couplings close to g t ∼ 2 could be excluded in the resonant region (m Y > 2m X ) with m Y < 100 GeV. For larger mediator masses, the cross section starts to fall due to the reduction of the phase space. In the off-shell region (m Y < 2m X ), the mono-Z cross section suffers from the same drastic decrease seen in figure 6 for the tt + / E T and monojet cases.
The same procedure can be repeated to constrain the parameter space of the model using mono-Higgs events on the basis of the results of the ATLAS search for dark matter production in association with a Higgs boson decaying into two bottom quarks [53]. This search results in a 95% CL upper limit on the visible cross section of 3.6 fb for a / E T threshold of 150 GeV. In order to estimate a limit, we generate events for (m Y , m X ) = (100, 10) GeV and g t = g X = 1, and require the two b-quarks to have a transverse momentum p b 1 T > 100 GeV and p b 2 T > 25 GeV, a pseudorapidity |η b | < 2.5 and to be separated in the transverse plane by an angular distance ∆R(b 1 , b 2 ) < 1.5. Moreover, we only select events exhibiting at least 150 GeV of missing transverse energy. We show again in figure 12 the / E T and leading b-quark transverse momentum distributions (blue lines) without and with applying the above-mentioned selection requirements. We then include a b-tagging efficiency of 60% and extract an upper limit on the g t coupling by comparing our results to the ATLAS limit. Coupling values of g t > 2 are found to be excluded for m Y > 2m X with m Y < 100 GeV. All other parameter space regions suffer from the same limitations as the mono-Z case. From our naive parton-level analysis, we have seen that mono-Z and mono-Higgs signals show promising signs of setting constraints on the parameter space of the model and therefore deserve dedicated studies, which will be reported elsewhere (see also ref. [71]). The sensitivity to such signals will benefit from applying more aggressive / E T thresholds to ensure the reduction of the corresponding backgrounds. As seen in figure 12, we obtain Resonance search constraints on top-philic DM Figure 13. Resonance search constraints from the LHC results at a collision centre-of-mass energy of 8 TeV on the simplified top-philic dark matter model presented in terms of the mediator mass m Y and the g t coupling. The different coloured areas are excluded by the diphoton [55] (orange), tt [56] (magenta) and tttt [57] (blue) searches. We include information on the mediator width to mass ratios (green curves). We assume a negligible branching ratio to the invisible sector. a rather hard / E T distribution [37], especially for mono-Z production. The result implies that an increase in the / E T threshold requirement in future analyses could lead to a significant improvement of the sensitivity, especially given the the fact that Standard Model backgrounds rapidly fall off with the increase in missing energy.

Dijet and diphoton resonances
Dijet and diphoton resonance search results could (in principle) be used to constrain the simplified top-philic dark matter model. Due to double-loop suppressions, mediatorinduced contributions to dijet and diphoton production are only relevant in the parameter space regions where m Y < 2m X , 2m t (i.e. where the mediator cannot decay into top quarks and/or dark matter particles). The partial mediator decay rate into gluons is then always dominant (as mentioned in section 2) since All LHC dijet resonance searches focus on the dijet high invariant-mass region, leading to no useful constraints on the top-philic dark matter model. The lowest mediator mass that is probed is ∼ 500 GeV, with a visible cross section restricted to be smaller than 10 pb [72].
Although the branching ratio of the mediator into a photon pair is very small, the background associated with a diphoton signal is low so that one expects to be able to obtain stringent constraints on the model from the diphoton search results. We focus here on the CMS 8 TeV diphoton search [55] that investigates resonance masses ranging from 150 GeV to 850 GeV and derives limits on the corresponding cross section. For instance, the 95% CL upper bound on the mediator-induced diphoton production cross section σ(pp → Y 0 → γγ) is of 20 fb (4 fb) for a mediator mass of 150 GeV (300 GeV). Making use of the pp → Y 0 cross section values shown in figure 6 and the Y 0 → γγ branching ratio computed from the formulas shown in section 2, we present diphoton constraints on the model in the (m Y , g t ) plane in figure 13. These results assume that the dark matter particle is much heavier than the mediator that can thus not resonantly decay invisibly. The constraints are found to be stringent below the 2m t threshold, where the g t coupling cannot be larger than 0.6.

Top-antitop resonances
For scenarios with mediator masses above the top-antitop threshold (m Y > 2m t ), tt resonance searches [56,73] can be used as probes of the model. In our setup, loop-induced resonant mediator contributions can indeed enhance the tt signal, in particular when there is a large coupling hierarchy (g t g X ) or mass hierarchy (2m t < m Y < 2m X ). We derive constraints on our model from the ATLAS 8 TeV tt resonance search [56] that relies on the reconstruction of the invariant mass of the top-quark pair to derive a 95% CL exclusion on the existence of a new scalar particle coupling to top quarks. The associated cross section limits range from 3.0 pb for a mass of 400 GeV to 0.03 pb for m Y = 2.5 TeV, assuming that the narrow width approximation is valid with a mediator width being of at most 3% of its mass and that there is no interference between the new physics and Standard Model contributions to the tt signal.
Constraints are computed using the NNLO mediator production cross section (see figure 6) and the relevant top-antitop mediator branching ratio derived from the formulas presented in section 2. The latter is in fact very close to one in the relevant region, the mediator decays into dark matter particle pairs being kinematically forbidden and those into gluons and photons loop-suppressed. The results are presented in the (m Y , g t ) plane in figure 13. This shows that scalar mediators with masses ranging from 400 GeV to 600 GeV could be excluded for g t couplings in the [1,4] range, the exact details depending on m Y and on the fact that the narrow-width approximation must be valid. This demonstrates the ability of the tt channel to probe a significant portion of the m Y > 2m t region of the model parameter space. In the region where 2m t , 2m X < m Y , the partial decay Y 0 → XX reduces the tt signal and therefore limits the sensitivity of the search.

Four-top signals
Scenarios featuring a mediator mass above twice the top-quark mass can be probed via a four-top signal, since the mediator can be produced in association with a pair of top quarks and further decay into a top-antitop system. Theoretically, the Standard Model four-top cross section has been calculated with high precision [74], but the sensitivity of the 8 TeV LHC run was too low to measure the cross section. Instead, an upper limit on the cross section at a centre-of-mass energy of 8 TeV has been derived [57,75]. The four-top production rate is constrained to be below 32 fb [57], a value that has to be compared to the Standard Model prediction of about 1.3 fb. Only models with new physics contributions well above the background (see e.g. ref. [76]) can therefore be constrained by the four-top experimental results.
In our top-philic dark matter model, the new physics contributions to the four-top cross section can be approximated by the ttY 0 cross section, the branching ratio B(Y 0 → tt) ∼ 1. Using the NLO cross section (see figure 6), we derive limits that we represent in the (m Y , g t ) plane in figure 13. A small region of the parameter space with g t > 2.5 and in which the mediator mass lies in the [2m t , ∼ 450 GeV] mass window turns out to be excluded. The weakness of the limit is related to the steeply decreasing cross section for pp → Y 0 tt with the increase in m Y .

The mediator width
In all the above studies where the final state does not contain any missing energy, the mediator width has been assumed narrow. Concerning the diphoton channel, this assumption holds within the entire excluded region as only loop-suppressed gluon and photon mediator decays are allowed. In the region where m Y > 2m t , the width of the mediator rises quickly with its mass, and the width over mass ratio rapidly exceeds the 3% value that has been imposed in the ATLAS tt resonance search [56] as can be seen in figure 13. The reinterpretation of the ATLAS results to a generic tt resonance model should therefore be made carefully, as the limit cannot be necessarily applied to scenarios featuring significantly larger mediator widths. This is shown in figure 13 by a dotted line, and we can also observe that most of the points that would have been excluded by the ATLAS search do not fulfil the requirement of a width below 3% of the mediator mass. In our excluded region of the parameter space, we allow the mediator width to reach 8% of its mass, by the virtue of the experimental resolution on the invariant mass of the tt system. This leads to the exclusion of scenarios with mediator masses up to 600 GeV.
The ATLAS resonance tt study claims that varying the width of the resonance from 10% to 40% for the massive gluon model results in a loss in sensitivity by a factor 2 for a 1 TeV resonance. An extension of the reinterpretation of the ATLAS limits on our simplified top-philic dark matter model to the case of larger resonance widths could then be performed by rescaling the limits by the appropriate correction factor. We have nonetheless found that no additional points are excluded even without rescaling the sensitivity of the search as the ATLAS analysis rapidly loses sensitivity for resonance masses above 600 GeV. Considering model points with a mediator width to mass ratio of at most about 8% therefore provides a realistic exclusion over the entire model parameter space.

Concluding remarks on direct mediator searches
Mediator resonance searches at 8 TeV show good prospects of constraining our simplified top-philic dark matter model, especially in the mediator mass range of 150-345 GeV and 400-600 GeV by means of the diphoton and top-pair searches respectively. So far, the tt resonance searches are strictly applicable to a limited parameter space region of the simplified model, and considering larger widths in the interpretation of the future results would allow for a more straightforward reinterpretation of the limits to a wider range of parameters. Concerning the four-top analysis, it can presently only exclude a restricted part of the parameter space, but future measurements are expected to lead to more competitive bounds.
Finally, the pp → tt + j channel could also be used to probe dark matter models coupling preferably to top quarks. This has been for instance shown in ref. [77] where a loop-induced production of ttj can in some cases lead to interesting constraints on topphilic models of new physics. In our case, they are nonetheless not expected to give more stringent constraints than the tt resonance searches. One could also consider the pp → tttj and pp → ttW t processes [77]. Because of the magnitude of the electroweak couplings, these processes are characterised by smaller cross sections than when four top quarks are involved, and are hence not likely to set more stringent constraints on the class of models under consideration.

Combined constraints
The final segment of our comprehensive study of top-philic dark matter simplified models is a combined study of astrophysical and collider constraints. We find that in the region where g X , g t ≤ π, the 8 TeV collider results that provide relevant bounds (once the relic density and direct detection constraints are imposed) originate from direct mediator production searches when the mediator further decays into a pair of Standard Model particles. Figure 14 illustrates our results and shows the scenarios that are excluded by resonant diphoton and top-pair searches as well as by the four-top analysis. All points in the plot accommodate the dark matter relic density and direct detection constraints, while the colours indicate points excluded by individual complementary collider bounds. The vast majority of excluded points lie in the region where 2m X > m Y with m Y ∈ [150, 600] GeV. This is the region where the mediator decay into a pair of dark matter particles is kinematically forbidden, ensuring large branching fractions for decays into Standard Model particles. The diphoton resonance search excludes points below the 2m t threshold, while tt results constrain the 400 < m Y < 600 GeV region. The four-top probe is able to exclude a narrow parameter space region close to m Y ∼ 2m t , in agreement with the findings shown in figure 13.
Relaxing the requirements on the relic density, the direct detection and the upper bound on the coupling strengths allows for another meaningful study of combined collider constraints. For this purpose we have performed a joint analysis of collider bounds on the top-philic simplified dark matter model in the scope of a four-dimensional parameter scan with a flat likelihood function over all dimensions. We have performed the scan by restricting the couplings to be smaller than 2π, as well as by allowing the mediator widths to reach 50% of the mediator mass. Figure 15 shows our results, where the upper left panel shows the model points excluded by the combination of all collider results, and the rest of the panels show the points excluded by individual LHC Run I collider results. We find that the 8 TeV monojet searches exclude model points which lie mainly in and around the triangle bounded by the m Y = 2m X and m Y = 2m t lines, where the characteristic g t which is excluded by the 8 TeV results is of O(10). The region in which the excluded LHC constraints on top-philic dark matter Figure 14. Results of our four-dimensional parameter scan projected onto the (m Y , m X ) plane once constraints set from the LHC results are imposed. The points excluded by the diphoton, the tt and the four-top considered searches all satisfy the relic density, narrow width and direct detection constraints.
points are located is reasonable, as we expect any significant monojet signal in the region where m Y > 2m X . Furthermore, we expect the branching ratio to missing energy to be lower in the region where m Y > 2m t due to the kinematically allowed decays into a pair of top quarks. This in turn leads to a lower signal cross section in all channels with missing energy and hence a lower number of points which can be excluded by monojet searches in the m Y > 2m t region.
The points excluded by the 8 TeV tt + / E T measurements lie in roughly the same region as the points excluded by the monojet search, but with a more defined edge of m Y = 2m t . Conversely, the 8 TeV tt resonance search provides constraints in the region of m Y ∈ [400, 600] GeV and m X 100 GeV, and is able to rule out g t couplings of O(1). The four top searches constrain roughly the same region of the (m Y , m X ) parameter space as the tt searches. However, the characteristic size of the couplings four top searches are able to constrain is significantly larger than the case of tt.
Finally the diphoton resonant search excludes m Y ∈ [150, 2m t ] GeV with 2m X > m Y , ruling out g t couplings larger than 0.6. In the (m Y , m X ) plane, we can observe that the constraints arising from all mediator resonance searches, i.e. the diphoton and tt analyses, are largely complementary to those issued from searches in channels with large missing energy. The results assume couplings smaller than 2π and Γ Y /m Y < 0.5, with no constraints from astrophysics or cosmology being imposed. In case of resonant tt, four top, and the combined constraints, we only show the 95% CL exclusion as the tt and four top results have not been obtained using a recast LHC analysis.

750 GeV diphoton excess
In the light of the possible new physics signal observed in the 13 TeV ATLAS [28] and CMS [29] diphoton data, we investigate whether the simplified top-philic dark matter model considered in this work can accommodate the features of the observed excess. The excess can be interpreted as a possible signal of a new particle with where σ γγ ≡ σ(pp → Y 0 → γγ). Near the resonance, the diphoton cross section is analytically approximated by [78] σ γγ (13 where C gg (13 TeV) = 2137 is the gluonic parton luminosity factor and Γ(Y 0 → gg) and Γ(Y 0 → γγ) are the partial decay widths given in eqs. (2.8) and (2.9). Eq. (6.2) should in principle also contain a contribution from the γγ initial state. Here we neglect the photon fusion production mechanism due to fact that in the top-philic dark matter model, the branching ratio to photons is always suppressed compared to the branching ratio to gluons (see eq. (4.2)). The relationship between the strength of the mediator couplings to gluons and photons is also one of the main differences between our simplified top-philic dark matter model and other dark matter models that have been proposed to explain the 750 GeV diphoton resonance excess. In the latter, the mediator couplings to gluons and photons are typically treated as independent parameters [79][80][81][82].
The different contributions to the mediator width when m Y = 750 GeV always include the tt, gg and γγ final states, while the partial decay into a pair of dark matter particles is subject to the value of m X . In order to determine whether our top-philic dark matter model can explain the diphoton excess, we hence only have to address two distinct regions of the parameter space.
• m Y < 2m X : In this region, Γ Y is obtained by summing the contributions of the decays into gg, tt and γγ final states. As the top decay channel is kinematically open, it will always dominate over the loop-suppressed gg and γγ modes, leading to Γ Y ≈ Γ(Y 0 → tt). The mediator-induced diphoton rate at a centre-of-mass energy of 13 TeV is then a function of the single parameter g t , where we have fixed m Y = 750 GeV. In order to reproduce the σ γγ 1 fb value that is necessary to explain the excess, one would hence naively need g t ∼ 19, which is way above the unitarity bound. Even without considering the mediator width, cosmology or astrophysics, we find that the top-philic dark matter model is not able to explain the diphoton excess when m Y < 2m X .
• m Y ≥ 2m X : The total mediator width in this region is well approximated by summing over the contributions originating from the decays into tt and XX pairs. The main difference compared to the previous case consists of the decay channel into a XX pair that is now kinematically allowed, which implies possibly suppressed branching ratios to the other final states. As the contributions of Y 0 invisible decays to a pair of dark matter particles only appear in the denominator of eq. (6.2), the maximum possible σ γγ cross section value is reached when Γ(Y 0 → XX) ∼ 0, i.e. when g X ≈ 0 or m X ≈ m Y /2. The resulting maximal diphoton rate then turns out to be identical to the one of eq. (6.3) so that the observed excess cannot be accomodated in our model.
The top-philic simplified dark matter model that we consider cannot accommodate the diphoton excess in any region of the model parameter space. Finally, we show the actual values of σ γγ (13 TeV) for the scenarios that feature relic density and direct detection cross section in agreement with data in figure 16, after restricting our selection to points featuring 730 GeV ≤ m Y ≤ 770 GeV. The largest cross section values that are found are at least two orders of magnitude too low in order to be able to accomodate the diphoton signal.

Conclusions
We presented a comprehensive analysis of simplified top-philic dark matter models, in the scope of collider physics, astrophysics and cosmology. Our study considered the full four dimensional model parameter space, where we treated the experimental constraints on the model space both separately and in conjunction with each other. The requirement of predicting the measured relic density Ω DM h 2 gives the most stringent constraint on the viable regions of the parameter space. Most of the region where m Y > m X cannot accommodate the observed relic density, except near the resonance m Y ∼ 2m X and for m X > m t . Direct detection data complementary excludes large portions of the parameter space in the m Y < m X region once experimental results from LUX and CDMSLite are accounted for. In the context of dark matter indirect detection, we studied prospects for further model constraints from gamma-ray flux measurements originating from dwarf spheroidal galaxies and the gamma-ray lines issued from the inner galactic region. In the specific model we consider, the dark matter annihilation cross section is p-wave suppressed, leading to indirect detection bounds which are too weak to provide additional constraints on the parameter space.
Collider searches from LHC Run 1 at √ s = 8 TeV can constrain the parameter space beyond the limits obtained from the relic density and direct detection, but apply mostly in the limit of coupling values 1. We found that for couplings of π, the resonant tt and diphoton searches are able to exclude a fraction of model points in the regions of m Y ∼ 400 − 600 GeV and m Y ∼ 150 − 350 GeV respectively, even upon assuming astrophysical and relic density constraints.
In addition to studying collider signatures of the top-philic dark matter simplified model as a complementary way of dark matter detection, we performed a study of collider constraints without assuming relic density and direct detection (as well as extended the parameter range to include coupling values of < 2π and Γ Y ≤ 0.5m Y ). Our results for a four dimensional parameter scan show that (in the scenario where astrophysical and cosmological constraints are not relevant), / E T + j and / E T + tt 8 TeV results provide meaningful bounds on the model parameter space in the 2m X < m Y < 2m t region, but only for g t , g X π. In the m X > m t region, the resonant tt searches are again able to exclude some model points in the m Y ∼ 400 − 500 GeV region, while γγ measurements provide constraints in the m Y < 2m t region. We have also explored the prospects of using rarer processes such as four-top production as well as mono-Z and mono-Higgs production to constrain our model. While we have not performed a detailed analysis we have found that these processes show promising signs of further constraining the parameter space of our model and deserve dedicated studies.
For the purposes of our study we have recast the CMS monojet and / E T + tt searches in the framework of MadAnalysis 5, which allows us to reliably extract constraints on our model, and can benefit future collider studies which go beyond our simplified model and even beyond dark matter searches. Another important aspect of our work, is the use of NLO QCD predictions for the / E T + tt process to constrain our model. While we find that K-factors for this process are close to one, the importance of taking higher order effects into account lies in the reduced theoretical uncertainties of the NLO results. We have shown that the uncertainties in the CL estimates significantly reduce with the inclusion of higher order QCD terms which clearly illustrates the importance of higher order corrections on the interpretation of dark matter searches at colliders.
In the context of the recently observed excess in the ATLAS and CMS measurements of the diphoton invariant mass spectrum, we consider the possibility that Y 0 decays to photons explain the excess of events around m γγ = 750 GeV. We find that due to top-loop suppressed couplings to gluons and photons, only non-perturbative values of g t suffice to fit the features of the excess. The work presented in this paper also represents a proof-of-concept for a unified numerical framework for dark matter studies at the interface of collider physics, astrophysics and cosmology in a generic model. As a part of consistency checks, we have ensured that the scan covers similar regions of the parameter space both in case of MadDM and micrOMEGAs. Fig. 19 shows the results for distributions of masses and couplings in the scans, where the blue/red lines refer to MadDM/micrOMEGAs respectively. Similarities in the distributions of fig. 19 indicate that parameter scanning was performed consistently between the two codes.

B.2 Details on the dark matter annihilation cross sections
In this Appendix we give the detailed analytic expression of the three annihilation processes described in Sec. 3.1. The s-channel annihilation cross section XX → tt (process (I)) is given by: Process (II) denotes the annihilation of dark matter into a pair of gluons via s-channel and is given by: Finally the process (III), namely XX → Y 0 Y 0 via t-channel is given by: where t 0,1 are the integration extrema: with t and u Mandelstam variables such that u = 2m 2 X + 2m 2 Y − s − t. In general the thermally averaged cross section can be approximated in the non relativistic regime by expanding the cross section in powers of the dark matter relative velocity v rel , with s m 2 X (4 − v 2 rel ), weighting with the appropriate Maxwell-Boltzmann distribution: where the index i indicated the annihilation process, x ≡ m X /T and T is the temperature of the dark matter gas. In case of s-channel annihilation, along the resonance the thermal average is much more complex and requires the full computation of the integral dx σv rel (x). The approximation given in eq. (B.6) holds in all regions far away from the resonance and is useful to show the dependence on v rel of σv rel (x) for each specific process.  For all processes (I), (II) and (III) the first coefficient is always null, A i = 0. The first non negligible term in the expansion eq. (B.6) is then B: This is equivalent to say that all three process are p-wave suppressed for dark matter annihilation at present epoch. The case of Dirac dark matter particles communicating with the SM via a pseudoscalar mediator has been described in [83], where analytic expressions for σv rel can be found. Similarly to scalar mediator Y 0 the t-channel process is again p-wave suppressed, while the s-channel annihilation is dominated by s-wave.

C Recasting of LHC searches within the MadAnalysis 5 framework
In this appendix, we detail the implementation, within the MadAnalysis 5 framework [20][21][22], of the two dark matter searches that we have investigated in this work. More precisely, this consists of the CMS-B2G-14-004 analysis [59] that probes final states comprised of a top-antitop system produced in association with a pair of invisible dark matter particles (see Section C.1) and the CMS-EXO-12-048 analysis [51] related to the production of a pair of dark matter particles together with a hard jet (see Section C.2). Both recasting codes have been validated within the version 1.3 of MadAnalysis 5, although the monojet search reimplementation is also compatible with the version 1.2 of the program. The simulation of the detector response is performed with the standard Delphes 3 package that we have run from the MadAnalysis 5 platform. In the monojet case, we have used the standard CMS detector parameterisation that is the shipped with MadAnalysis 5, while in the top-antitop plus missing energy case, we have designed a dedicated detector card. For both setups, jets are reconstructed on the basis of the anti-k T algorithm [84] with a radius parameter set to 0.5, as implemented in FastJet [85].
The validation of both our reimplementations is based on material provided by CMS. Two UFO models [86], one for each of the recast analyses, have been shared so that we have been allowed to generate specific dark matter signals for which CMS has released public cutflow charts and differential distributions. Using MG5 aMC [17] (with the leading order set of CTEQ6 parton densities [87]) and Pythia 6 [88] (with the Z * 2 tune [89] for the description of the underlying events) for the simulation of the hard scattering process and of the parton showering and hadronisation, respectively, we have generated signal events that have been analyzed with MadAnalysis 5. Our results have been confronted to the CMS official numbers, which has allowed us to assess the validity of our recasting codes. Our simulation procedure moreover includes the generation of matrix elements containing up to two extra jets that we have merged according to the MLM prescription [90,91], the merging scale being set to 40 GeV.

C.1 The CMS top-antitop plus missing energy CMS-B2G-14-004 search
In order to validate our reimplementation of the CMS-B2G-14-004 search in MadAnalysis 5, we focus on a new physics model that features the production a pair of dark matter particle X of mass m X = 1 GeV in association with a top-antitop pair via a four-fermion interaction. The CMS event selection strategy requires a large amount of missing transverse energy, a single isolated lepton and multiple jets, and uses 19.7 fb −1 of proton-proton collision data recorded at a center-of-mass energy of √ s = 8 TeV.
The CMS-B2G-14-004 analysis relies on single electron and muon triggers, with lower p T thresholds of 27 GeV and 24 GeV respectively, and the reconstructed electron (muon) candidate is imposed to be isolated in such a way that the sum of the transverse momenta of all objects lying in a cone of radius R = 0.3 centered on the lepton has to be smaller than 10% (12%) of the lepton p T . Event preselection finally requires that the lepton p T is larger than 30 GeV and pseudorapidity |η| is smaller than 2.5 (2.1 for muons). It additionally demands the presence of at least three jets of p T > 30 GeV and |η| < 2.4 with one of them being b-tagged, as well as missing energy / E T > 160 GeV. The signal region is defined by selecting events with a large amount of missing transverse energy / E T > 320 for which the transverse mass M T that is constructed from the lepton and the missing energy is larger than 160 GeV. Moreover, the missing transverse momentum and the two leading jets are asked to be well separated in azimuth, ∆Φ j 1,2 , / E T > 1.2, and the M W T 2 variable [94] is enforced to be greater than 200 GeV.
In Table 7, we confront the cutflow chart that has been obtained with MadAnalysis 5 to the official results of CMS for the benchmark scenario under consideration. For each where n i and n i−1 mean the event number after and before the considered cut, respectively. The relative difference information given in the table corresponds to the difference between the MadAnalysis 5 and the CMS efficiencies, normalized to the CMS result, -43 -An agreement at the percent level has been found all over the selection procedure. Moreover, we compare several (normalized) differential distributions as calculated with Mad-Analysis 5 when all selection steps but the one related to the represented kinematic variable are included with the public CMS results in figure 22. A very good agreement can again be observed.

C.2 The CMS monojet CMS-EXO-12-048 search
The validation of our implementation of the CMS-EXO-12-048 search in MadAnalysis 5 has been achieved on the basis of a benchmark scenario that is inspired by Refs. [95][96][97][98]. In this context, monojet events arise from the associated production of a pair of invisible Dirac fermions of mass of 1 GeV with at least one hard jet. The interactions of the dark particle with the Standard Model are mediated by a new gauge boson Z of mass and width of 40 TeV and 10 GeV respectively, and all new physics interactions have been assumed to have a vector coupling structure and a strength set equal to 1. Concerning our signal simulation setup, we have imposed that all parton-level jets have a transverse momentum p T larger than 20 GeV and that the leading jet has a p T > 80 GeV.
The CMS monojet search relies on an integrated luminosity of 19.7 fb −1 of protonproton collisions at a center-of-mass energy of √ s = 8 TeV. It focuses on a signal containing a very hard jet with a transverse momentum satisfying p T > 110 GeV and a pseudorapidity smaller than 4.5 in absolute value. A second jet is moreover allowed, provided that its transverse momentum is larger than 30 GeV, its pseudorapidity satisfies |η| < 4.5 and if it is well separated from the first jet by 2.5 radians in azimuth. Events featuring more than two jets (with p T > 30 GeV and |η| < 4.5), isolated electrons or muons with a transverse momentum p T > 10 GeV or hadronically decaying tau leptons with a transverse momentum p T > 20 GeV and a pseudorapidity satisfying |η| < 2.3 are discarded. The analysis then contains seven inclusive signal regions in which the missing energy / E T is required to be above specific thresholds of 250, 300, 350, 400, 450, 500 and 550 GeV respectively.
The selection strategy of the CMS monojet analysis thus consists of six preselection cuts followed by one region-dependent cut, when we ignore the first two requirements of the analysis related to the cleaning of the events from the detector noise that cannot be handled with Delphes 3. For the benchmark scenario under consideration, we compare the results that have been derived with our MadAnalysis 5 reimplementation with those provided by the CMS collaboration in Table 8.
We have found that all selection steps are properly described by our implementation, with the exception the missing energy selection / E T > 250 GeV for which a disagreement of about 20% has been observed. It is however not uncommon that low missing energy is difficult to simulate with a fast-simulation of the detector based on Delphes 3. We have verified that for missing energy values of interest, the description of the missing energy agree relatively well with CMS, as illustrated in figure 23 where we compare, for a benchmark scenario where the Z mass has been set to 900 GeV, the missing energy distribution as obtained by CMS to the one derived with MadAnalysis 5.