Contact interactions and top-philic scalar dark matter

We investigate the phenomenology of a scalar top-philic dark matter candidate when adding a dimension-five contact interaction term, as motivated by possible underlying extensions of the Standard Model such as composite Higgs models. We show that the presence of contact interactions can have a major impact on the dark matter relic density as well as on its direct and indirect detection prospects, while the collider phenomenology of the model is unaffected. This underlines the complementarity of collider and cosmological constraints on dark matter models.


Introduction
That about 27% of the universe's energy budget is made up of matter which cannot be described by the Standard Model (SM) is one of the foremost mysteries in particle physics. Dark matter (DM), so named for being non-luminous and non-absorbing [1], is non-relativistic matter which does not behave like the baryonic matter of the SM. The nature and origin of DM is still unknown, despite compelling evidence [2] for its existence, and numerous physics programs have been established in an attempt to detect it. These include direct and indirect detection attempts and efforts at colliders such as the Large Hadron Collider (LHC) at CERN, none of which have yet made conclusive discoveries.
The "cold" DM theory postulates that DM has been (and continues to be) nonrelativistic since the beginnings of galaxy formation, with its evolution governed by the Boltzmann equation. In the early universe, DM was in thermal equilibrium with SM particles. As the universe expanded and thereby cooled, the DM collision rate dropped to the Hubble expansion rate, at which point the particles "froze out" and decoupled. The observed DM density, or relic density Ω DM h 2 , became constant. Candidates for cold non-baryonic DM are therefore strongly constrained by the relic density, which has been measured to be [3] Ω DM h 2 = 0.1186 ± 0.0020, (1.1) and which is controlled by the annihilation cross section of the DM candidate. A number of candidates for DM have been proposed, but the WIMP paradigm, in which a new dark state couples to the SM through a generic weak interaction and has a mass ranging from several GeV to a few TeV [2,4], persists as a promising avenue.
In this work we build on existing models of a scalar heavy DM candidate S that manifests as a weakly interacting massive particle (WIMP) coupling to the SM top quark t via a Yukawa-type term which involves a heavy fermionic mediator T . Moreover, we constrain the mass m S of our real scalar DM particle S to be larger than that of the top quark (m t ), focusing on candidates S with 200 GeV m S 3 TeV. As we are keeping in mind that S and T could be resonances emerging from a composite Higgs theory, we require their masses to lie within an order of magnitude of each other. In this m S range the DM annihilation channel SS → tt dominates, and with m S well separated from the m t threshold effects which would arise in more mass-degenerate regimes can be avoided.
The DM resonance here considered is thus top-philic, as suggested in composite extensions of the SM in which the top quark plays a special role, and couples to the SM via a fermionic mediator T which is also heavy. This occurs via a t-channel interaction originating from an ST t operator. Many scalar DM models including a t-channel fermionic mediator have previously been studied [5][6][7][8][9][10][11], as top-philic DM models [12][13][14][15][16][17][18] and heavy DM with masses ranging up to several TeV have been discussed in the literature [19]. In particular, recent investigations [9,10] have shown that while leading order (LO) calculations allow for a heavy DM, next-to-leading-order (NLO) corrections to the annihilation cross section lead to significant modifications of the existing constraints on the model parameter space.
In the current study, we envision that the DM candidate S and the fermionic mediator T may arise within a composite Higgs model with an underlying fermionic construction as composite bound states. Such a composite bound state would be expected to have a mass of the order of the energy scale of the theory Λ, which could be expected to be several TeV [20]. While the possibility that DM arises from a composite Higgs model as an additional pseudo-Nambu-Goldstone boson (pNGB) has been well studied, the prospect that it is a heavy resonance has received less attention. It is this structure which we bear in mind throughout this work, although the couplings and parameters are left free to allow for the discussion of a more general situation. With this aim of generality, we recall that dimension-five operators are a generic feature of a broad range of Beyond the Standard Model (BSM) theories. In particular, these include composite Higgs models arising from strong dynamics, where higher dimensional operators do not decouple [20] and may therefore be relevant at colliders and in direct and indirect detection experiments. Instead of focusing on a particular theory, our methodology complements and generalises earlier studies relying a simplified model construction [9,10] and add to such a modelling an additional dimension-five operator SStt which emerges from a contact term between the S and t states. This independent dimension-five operator comes with an unknown O(1) Wilson coefficient, whose sign and magnitude could substantially modify the existing limits obtained through only dimension-four operators.
In fact, the addition of this term contributes to the relic density calculations by opening an area of parameter space not previously allowed, where the Yukawa termỹ t governing the ST t operator is very small. The interplay between the dimension-five operator and Yukawa term is fully parametrisable according to the masses of the particles, but is subject to possibly large interference for m S ≈ m T . We begin this work by investigating the predictability of the behaviour of the system across mass compression scales and with the inclusion of the additional dimension-five term. We have determined the analytical function which predicts the interplay between the dimension-five operator and the Yukawa coefficient which results in the correct relic density, potentially bypassing the need for intensive numerical calculations in order to determine the parameters yielding the correct relic density. In doing so we also display the need to account for co-annihilations in the highly compressed regime.
We then calculate expected bounds from DM direct detection experiments from the DM candidate S scattering off atomic nuclei. While the absence of a valence top quark from the nucleus makes the DM-gluon interactions the only avenue of detection, we may expect some suppression due to the loop-generated nature of these interactions. We see in fact that the addition of the new interaction term greatly improves the direct detection prospects, where many potential model configurations live above the neutrino floor and are therefore hypothetically accessible. Additionally, some models are within reach of the XENON-1T experiment. In these discussions we consider both possible signs of the Wilson coefficient. In studying the indirect detection constraints and prospects, we also find that more of the parameter space becomes accessible. Finally, we examine the collider constraints relevant to our model, considering the reinterpretation of existing mono-jet and multi-jet analyses to assess our model against current bounds. We find no improvement on existing collider constraints due to the dimension-five effects, and instead update the bounds given in previous analyses in the light of full LHC run 2 results. We moreover estimate the sensitivity of the future high-luminosity phase of the LHC to the considered class of top-philic scalar DM models.
In the following we begin with a description of the theory in section 2, defining the model and cross sections relevant to the relic density calculation. In order to map out the parameter space of interest, the relic density is calculated and relevant values of the Yukawa and dimension-five couplings are determined. In section 3 we study the impact on the direct detection constraints by the addition of this dimension-five term, while indirect detection implications are discussed in section 4. In section 5 we describe the hadron collider phenomenology of the model. We summarise our work in section 6 and we detail in appendix A the fit of the relic density that we have implemented. This exhibits the interplay between the variables of the model to a parametrisable curve, allowing for the prediction of relevant couplings according to the mass of the DM and mediator while accounting for co-annihilations.

Heavy dark matter with a t-channel mediator
We consider a simplified model description of the top-philic DM model [6,7,9,10] in which the SM singlet scalar DM candidate S couples to the SM through a contact interaction with the SM Higgs doublet φ, and a Yukawa-type interaction with the top quark and a vector-like top partner T . This vector-like partner has the same SM quantum numbers as the right-handed top quark, and we enforce that the masses of the new states satisfy m S ≤ m T . The model has a discrete Z 2 symmetry, under which S and T are odd while all SM particles are even, thus guaranteeing the stability of S. Following Refs. [9,10] In our notation,ỹ t stands for the ST t Yukawa coupling strength and λ for the strength of the Higgs portal to the dark sector. The final term in this Lagrangian is a dimension-five operator linking the DM particle S to the top sector via an effective contact interaction. The unknown coefficient of the contact term operator cannot be determined within the effective theory, and we parametrise it as C/Λ on dimensional grounds, where C is a dimensionless coefficient and Λ parametrises the scale at which the model is embedded into a more fundamental theory. 1 In order to solely focus on the top-philic nature of the model, we assume that the DM coupling λ to the Higgs field vanishes, departures from this hypothesis being discussed in Ref. [6]. Thus, the effective Lagrangian (2.1) yields a four-dimensional parameter space with two masses, m S and m T , and two couplings,ỹ t and C/Λ. In this work we will parametrise the mass-splitting between S and T by introducing the dimensionless quantity r = m T /m S − 1.
The DM and collider phenomenology of this top-philic DM model but without additional contact interaction has been discussed in Refs. [6,7,9]. The study in Ref. [9] demonstrated the importance of QCD radiative corrections, which have a major impact on the parameter space which reproduces the observed DM relic density, as well as on direct and indirect detection prospects. This shows that the DM phenomenology of the model is potentially sensitive to a priori suppressed corrections, and motivates studying the impact of the contact interaction parametrised by C/Λ.
Diagrams relevant for the calculation of the relic density through the SS → tt channel are shown in figure 1, where the first diagram contributes to the leading-order (LO) cross section and the next three to its next-to-leading-order (NLO) corrections in the strong coupling α s . To these diagrams we add the process which emerges as a result of the dimension-five contact term (last diagram in the figure). Following Ref. [10], the NLO annihilation cross section is well approximated by with r = m T /m S − 1, N c = 3 and C F = 4/3. The SStt contact interaction yields an additional contribution to the thermally averaged cross section whose leading component reads In principle, the relic density contains dimension-four contributions, dimension-five contributions, and their interference. In the relevant regions of the parameter space, the interference is however found to account for at most 1-2%. It is therefore neglected in the following. We determine the relic density numerically (and use these results for the DM and collider phenomenology analysis in sections 3, 4 and 5) and present the results first, before interpreting them through a semi-analytic description which is obtained by a fit to the numerically obtained solutions. The relic density is calculated by including both the NLO and SStt contributions to the annihilation cross section. This allows for the determination of the regions of the parameter space (m S , m T ,ỹ t , C/Λ) in which the relic of the DM candidate matches experimental data. In order to estimate the DM relic density including not only the NLO effects but also all relevant annihilation and co-annihilation channels, we employ the MicrOMEGAs [21] framework, for which we generate a CalcHEP [22] model file through its interface to FeynRules [23,24]. For the parameter scan we vary the mass m S (m T ) between 200 GeV and 3000 GeV (3500 GeV), and allow the Yukawa couplingỹ t to lie within the range [10 −4 , 6]. For the C/Λ coupling we then scan over the interval [10 −3 , 10 −5 ] GeV −1 .
The relative impact of the LO, NLO and contact operators to the cross sections is illustrated in figure 2, for a small sample of benchmark scenarios. We have considered three DM masses of m S = 500, 1000 and 2000 GeV and fixed the heavy-top mass m T through r = 0.6. We have then varied the Yukawa couplingỹ t from 0 to 4, and updated the C/Λ value to match the observed relic density. The figure exhibits the interplay between the different terms that contribute to the total annihilation cross section σv = σv N LO + σv SStt , where the separate contributions sum to the thermally averaged cross section required to yield the observed DM relic density. As can be seen, smallỹ t yields a small NLO cross section which is compensated for by the contact operator contribution. Also shown in the figure is the LO contribution, which is accounted for in the NLO result (2.2), where the VIB processes become increasingly relevant for higher masses m S . This figure shows the importance of the dimension-five contribution from the contact operator to the total annihilation cross section in the regime of small Yukawa coupling, and motivates further investigation into potential modifications to the phenomenology and direct and indirect detection constraints. The strength of the contact interaction is of order C/Λ 0.2 TeV −1 (or less, if the cross section is dominated by the NLO contribution) and in the following we determine a semi-analytic approximation for the value of C/Λ required to reproduce the observed DM relic density.
The DM relic density obeys the Boltzmann equation 2 where H is the Hubble constant, n denotes the number-density of S, n eq is the equilibrium number density, and σv is the thermally averaged annihilation cross section. An approximate solution to the Boltzmann equation is given by [25] Ω DM h 2 ≈ 1.04 × 10 9 M P l where M P l is the Planck mass, g * is the total number of effectively massless degrees of freedom, and x F is the ratio of the dark matter mass to the freeze-out temperature. The a and b coefficients are functions of the masses and are defined from the non-relativistic annihilation cross section, which we can expand as For the currently considered model, the thermally averaged cross section is a sum of equations (2.2) and (2.5), The dependence on the couplingsỹ t and C/Λ can be factored out, so that with A(m S ) being a function of the DM mass and B(m S , m T ) being a function of both the DM and the mediator masses. Therefore, using that Ω DM h 2 ∼ 1/ σv [26], solving for (2.12) The b (x F , g * (x F )) function is then determined from a fit to the numerical result. For the result in eq. (2.11) we neglected co-annihilation effects. They play a role when the DM candidate S and the mediator T are nearly mass degenerate, which is the case in parts of the parameter space considered. The co-annihilation effects, too, can be treated semi-analytically, and we refer the interested reader to Appendix A. Here, we just quote the final result that generalises eq. (2.11), where A(M s ) and B(m S , m T ) are given in eq. (2.12), r = m S /m T − 1, and the coefficient b (x F , g * (x F )) is determined by a fit to the numerical results to be b = 6.0 × 10 −9 ± 0.2 × 10 −9 GeV −2 . The other coefficients which parametrise the co-annihilation effects are fitted to (α, β, γ) = (0.4, 1.9 × 10 −4 , 6.2 × 10 8 ) for m S ≤ 1.2 TeV, and (α, β, γ) = (0.7, 3.0 × 10 −3 , 1.8 × 10 4 ) for m S > 1.2 TeV. The parameter γ has been raised to a dimensionless ratio featuring Λ = 3.5 TeV, which is the maximum value for m T used in the scan. Such a value hence provides an indicative scale of the effective model (such as the limit of validity of the theory or the scale of compositeness).

Direct detection prospects and bounds
As the DM candidate in our model interacts with ordinary matter, its properties can be probed at direct DM detection experiments. In such experiments, the collision of S with a nucleus of the detector material can leave hints to be observed through the recoil energy of the nucleus. The rate at which this occurs is related to the nucleon-DM cross section as predicted in our model. As the DM scalar S only interacts with the Higgs boson (although we neglect such interactions) and the top quark, it does not have any tree-level interactions with valence quarks of the nuclei, and DM-gluon scattering at the loop level is the dominant contribution to DM-nucleus scattering.
In evaluating the scattering cross section with a nucleon we follow Ref. [27], matching the effective theory describing the DM-nucleon interactions to the full theory through higher-dimensional operators. In the case of scalar DM, only the spin-independent cross section is relevant. As the low velocity of the DM leads to a small momentum transfer [28], the interaction between the DM and nucleons can be described by the following effective Lagrangian where the effective operator is defined at the mass scale of the mediator. The Wilsoncoefficient C g S receives contributions f d4 resulting from the renormalisable part of the Lagrangian (2.1), as well as contributions f d5 originating from the dimension-five contact interaction examined in this work, The full expression for f d4 (m S , m t , m T ) has been determined in Refs. [9,27] (and are in particular collected in the appendix of Ref. [27]). The new contribution f d5 arises from the addition of the SStt coupling to the model Lagrangian. It leads to an additional diagram to the total amplitude for the DM interaction with the nucleus, displayed in figure 3. The effective coupling is written as where "| GG " indicates that the coefficient for the terms proportional to G a µν G aµν have been extracted from the quark propagator in the gluon background in the limit of zero gluon momentum [29] Figure 3. Additional diagram to be considered in the evaluation of the S coupling to a nucleus via the constituting gluons of the latter, once the model Lagrangian includes an SStt dimension-five operator.
With the Wilson coefficient determined, the spin-independent DM-nuclear cross section σ A (for n p protons and n n neutrons in a nucleus of mass m A ) and DM-proton cross section σ p read For a nucleon N we moreover have with the quark mass fractions f Tq being given in table 1 of Ref. [27]. In this notation, the dependence on the strong coupling constant is implicitly absorbed in the definition of f N .
The scattering cross sections in eq. (3.5) are proportional to C g S 2 , and C g S receives a positive contribution proportional toỹ 2 t from the renormalisable interactions in the TeVscale theory. In contrast, the contact interaction contribution is proportional to −C/Λ. The sign of C is not fixed at the EFT level, such that the contact interaction can either increase or decrease the direct detection cross section, depending on the sign of C/Λ. While the sign of C/Λ plays an important role here, we recall that it was not relevant in the determination of the DM relic density.
In figure 4 we present the DM-proton cross section σ p across the full considered m S mass range. Recalling that we consider the parameter ranges C/Λ ∈ [10 −5 , 10 −3 ] GeV −1 andỹ t ∈ [10 −4 , 6], we expect non-zero contributions from both the NLO and dimensionfive processes in determining the full cross section displayed in the figure. Here and for the remainder of the paper, we consider r = m T /m S − 1 ∈ [0. 1,9], where the upper limit ensures that S and T are separated in mass by at most an order of magnitude. This is motivated by the idea that both emerge as resonances in a composite Higgs model. The lower limit is chosen as the compression scale at which co-annihilation effects start to become apparent (see appendix A for further details    1,9] interval and is derived from the value of the other model parameters and the relic density. The bottom limit of 0.1 ensures that we are in a well-understood regime of spectrum compression, and the upper limit is chosen such that the vector-like resonance and the scalar DM are not too different in mass. Additionally, we display the cross section obtained through use of the maximum possible Wilson coefficient |C|, which is displayed in dark blue. In those scenarios the dependence of the relic density onỹ t is negligible with respect to the contribution from the SStt interactions. Finally, the red dashed line represents the neutrino floor [30], the orange dashed line indicates the XENON 1T reach [31], and the red solid line shows the 90% confidence exclusion of the XENON 1T experiment [32].
certainly be modelled by the tools at hand, we choose to stay far from co-annihilation effects as we do not fully account for them.
At each DM mass, we first single out the ensemble of scanned scenarios featuring the right relic density, and then plot the σ p values associated with each scenario from a selected subset. In this subset of models, the absolute value of the Wilson coefficient |C/Λ| is minimum, and the value of the Yukawa couplingỹ t is thus maximum. The corresponding value for the compression factor r is also depicted, through the red-and-blue colour scheme. Moreover, we superimpose to our predictions constraints and the projected reach from the XENON experiment, as well as the neutrino floor. Conversely, the dark blue line includes the results obtained when selecting a subset of scenarios for which the Wilson coefficient is maximum, and the Yukawa couplingỹ t is thus at a minimum. This figure demonstrates that many viable models from the relic density standpoint lie above the neutrino floor, and that some models, particularly in the low DM mass region, are within the reach of future upgrades of the XENON experiment. Those are thus in principle testable in the future. In the figure, the difference in behaviour for positive and negative C/Λ is also highlighted.
For a positive dimension-five coupling, the viable regions of parameter space display significant overlap across r = m T /m S −1 values. This is in contrast to the case of a negative dimension-five coupling, where the resulting parameter space is more spread out, and many scenarios with lower r values result in very small σ SI p , which are not viable for detection.
The benchmark points with r 1 values are common to both setups, i.e. the reddish band across the centre of the figure is unchanged for positive and negative C. We recall that such a band corresponds to low |C/Λ| and largeỹ t , where C is non-zero and the Yukawa coupling has taken over entirely in the calculation of the relic density. The implications for direct detection are markedly different from the relic density calculations for low |C/Λ|; in the case of the relic density, the minimum value for |C/Λ| was chosen to reflect the case where the dimension-five term no longer modifies the relic density produced by the Yukawa term, producing a relic density comparable to the case C = 0. However, the same minimum value for |C/Λ| clearly modifies the direct detection prospects in a non-negligible way. As soon as |C/Λ| increases from its minimum value, so does the DM-proton scattering cross section to finally reach the bright blue line on the figure, that corresponds to a maximum |C/Λ| value. Here, the dimension-five addition completely dominates the DM-proton cross section and takes over from the Yukawa term, leading to an unchanged cross section under the C sign flip. Even for 'small' C/Λ, as is considered here, the contribution from the dimension-five operator is still larger than the one proportional to theỹ t Yukawa coupling. It is clear from the figure that the difference between positive and negative C emerges for smaller values of r. For r values smaller than about 1 (in particular for m S > 500 GeV), the C/Λ contribution is about the same order of magnitude as theỹ t contribution, such that when C/Λ is positive and the contributions have the same sign, they add constructively and make a larger cross section. When C/Λ is negative, they add destructively and σ p may be pushed down. We find that, for a representative point r = 0.1, heavy and light DM behave very differently; for light scalars the addition of the dimension-five greatly increases the cross section, but for positive and negative C there is no change. In the case of heavy scalars and r = 0.1, a sign flip in C leads to destructive interference. This contrasts with r ≥ 1 configurations in which the cross section is unchanged under the flip of the sign, no matter m S .
The limiting values for the intersection with the neutrino floor, the XENON 1T reach, and the 90% confidence exclusion of the XENON 1T, are given in table 1. In this table we give the mass m S at which the relevant exclusion intersects the distribution of phenomenologically viable scenarios, for both positive and negative Wilson coefficients. It is notable that the most significant difference arises in the case of the neutrino floor. There is a much larger number of scenarios below it, and thus potentially hard to probe, in the case of a negative Wilson coefficient. We now investigate this last property deeper.
In figure 5, we determine the region of interest in the (C/Λ,ỹ t ) parameter space for a number of illustrative benchmarks in increasing m S . The shaded regions correspond to parameters where the DM-proton cross section crosses the neutrino floor, and is therefore potentially reachable by experiments. The shaded region has been coloured light blue if the region leads to a relic density where the DM is over-abundant, and darker blue if the DM is under-abundant. The black line indicates the parameters yielding the correct relic density. Given that Ω DM h 2 ∝ σv −1 , the area outside/above the relic function line is under-abundant and not excluded, and the area below the relic line is over-abundant and therefore excluded. The dotted black line indicates scenarios with a correct relic density, but in a region where the DM-proton cross section lies below the neutrino floor. The C/Λ parameter clearly plays a dominant role in pushing the DM-proton cross section across the neutrino floor in the region corresponding to the correct relic density when only relying on these two terms. We can relate the bright blue line visible in figure 4 that corresponds the maximum C/Λ coupling to the horizontal portion of the black lines in figure 5. Moreover, the top plots in figure 5, both at m S around 200 GeV, illustrate how a change in mass splitting modifies the Yukawa coupling impacting the features of the models in direct detection experiments. A stronger compression indeed reduces the Yukawa coupling value at which the models will become visible relative to the neutrino floor.

Indirect detection prospects and bounds
Experiments which seek to detect DM through indirect methods aim to observe hints from DM annihilations or decays into SM particles that then reach us in the form of gamma or cosmic rays which may travel through the universe with little other interaction [33]. The flux of these SM particles depends on the annihilation cross section of the DM, the relative branching ratios of the different final-state particles produced in those annihilations, the mass m S , as well as astrophysical constraints. The indirect detection of DM through cosmic rays holds the advantage that we can potentially detect DM on galactic or cosmological scales. In order to assess whether the indirect detection bounds may differ from those in previous works [9], we again utilise FeynRules for the generation of UFO model files [34], this time as inputs into MadGraph5 aMC@NLO (MG5 aMC) [35]. Next, we simulate DM annihilation at close to zero velocity, using Pythia 8 [36] to describe parton showering and hadronisation.
In order to compare our predictions to the exclusions from experiments, we examine the gamma-ray spectrum from bb and tt DM annihilation in three cases of interest: with no dimension-five operator present, with the lowest allowed value for the dimension-five coupling in the range deduced from the relic density scan, C/Λ ∈ [10 −5 , 10 −3 ] GeV −1 , and the corresponding highest Yukawa coupling for given mass setups, and finally with the highest dimension-five coupling and lowest Yukawa for the given mass setups.   state, We hence use this relation to rescale Fermi-LAT expectations which are associated with the bb channel from dwarf spheroidal galaxy future data, when assuming 15 years of Fermi-LAT operation [38]. The second conclusion from figure 6 is that the dimension-five term in the model Lagrangian does modify the associated gamma-ray spectrum, and is therefore expected to impact on indirect detection constraints. In order to assess this further, we present in figure 6 (right) the σv cross section at zero velocity, indicating the contribution of the dimension-five operator to be added to the full NLO results presented in Ref. [9]. 3 We additionally superimpose to our results the constraints that could be extracted from Fermi-LAT and cosmic ray data and expectations. Concerning the latter, DM annihilations into tt systems can indeed be constrained with proton anti-proton cosmic ray data [37], given that the tt and bb spectra display the same behaviour as in figure 6.
The bounds achievable through the inclusion of both the dimension-five operator and NLO QCD contributions are shown in the summary plot of figure 9, discussed in our conclusions. We recall that annihilations into pairs of gluons are neglected, as they are only relevant below threshold, for m S < m t (i.e. a region into which we do not venture).

Collider phenomenology
Finally, and in addition to modifying the astrophysical constraints, the additional vertex due to the dimension-five operator may modify the collider constraints on the model. Experimental searches for DM form an important part of the new physics search program at the LHC, and there exist a number of previous physics searches which may be reinterpreted to constrain the model examined in this article. The potential collider signatures of the model include a mono-jet channel and a multi-jet plus missing transverse energy ( / E T ) channel, as well as a tt + / E T mode. Figure 7 shows examples of diagrams contributing to these signatures. For the mono-jet final state, diagram (1a) of figure 7 is independent of the dimension-five operator, while diagram (1b) contains one ttSS vertex. The amplitude of this latter diagram is thus proportional to C/Λ. Analogously, for the pp → ttSS process (that leads to a multi-jet plus missing energy final state once top decays are accounted for), diagrams (2a) and (2b) of figure 7 do not depend on the dimension-five operator, while diagram (2c) contains one ttSS vertex. The full pp → jSS and pp → ttSS cross sections can thus be expanded as where σ 0 stand for ('bare') cross sections in the top-philic DM model without the added dimension-five operator,σ int jSS is the contribution from the interference of diagram (1a) with diagram (1b),σ int ttSS is the contribution from the interference of diagram (2a) and (2b) with diagram (2c), andσ dim5 jSS andσ dim5 ttSS are the 'bare' cross sections arising solely from the amplitude described by diagrams (1b) and (2c). For all contributions we factored out the BSM couplingsỹ t and C/Λ, and indicated the dependence on the BSM particle masses m T and m S .
Obtaining the correct relic density imposes a bound of C/Λ 0.175 TeV −1 (see our conclusions and figure 10). The complementarity with cosmology therefore suppresses the dimension-five contributions in eq. (5.1). Given C ∼ O(1), we obtain effective scales of 5 − 100 TeV for the relevant parameter space regions which yield the correct relic density. As currently pursued searches for dark matter at the LHC typically probe scales of the order of 1 TeV or less, our predictions can safely be trusted in terms of the validity of the effective field theory. Moreover, jSS production yields negligible cross sections once a cut as used in mono-jet searches is imposed on the jet transverse momentum (p T ), both for the dimension-four and dimension-five components of the cross section. We are thus left to consider the ttSS final state. For a centre-of-mass energy of √ s = 13 TeV, the dimensionfive operator contributes to the cross section for at most 0.003 fb for m S > 200 GeV, so that the inclusion of the dimension-five operator does not impact the top-philic DM model's collider phenomenology. The latter instead relies on usual vector-like top partner production, as induced by QCD interactions and whose pair-production cross section lies deep in the fb regime probed at the LHC [39].
To assess the LHC constraints on the model, we generate an NLO UFO model with FeynRules [23,24,34,40], and then use MadGraph5 aMC@NLO [35] in conjunction with Pythia 8 [36] to generate hadron-level events for the considered process pp → TT → ttSS. In our simulation chain, T decays are handled with MadSpin [41,42] and the matrix elements are convoluted with the LO set of NNPDF 3.0 set of parton densities [43,44]. The simulation of the response of the LHC detectors and event reconstruction are performed with Delphes 3 [45] (with appropriate detector parametrisations), that internally relies on FastJet [46] and its implementation of the anti-k T algorithm [47]. In our framework both detector simulation and event reconstruction are dealt with MadAnalysis 5 [48][49][50][51], which is then used for the extraction of the CL s exclusions relative to recent ATLAS and CMS searches for dark matter in the multi-jet + / E T and in the tt + / E T modes. We evaluate multi-jet plus missing energy constraints by recasting the results of the ATLAS CONF 2019 040 [52] search that covers a luminosity of 139 fb −1 and targets final states featuring at least two hard jets in association with missing momentum. For the tt + / E T limits, we consider the CMS-SUS-17-001 search [53] that analyses 35.9 fb −1 of data and focuses on stop pair production in a final state comprising two opposite-sign isolated leptons, two hard jets, and well separated missing transverse energy. Both these analyses are available within the Public Analysis Database of MadAnalysis 5 [54,55]. Although the CMS-SUS-17-001 analysis has recently been updated to 137 fb −1 of data [56], such an update is not yet included in the MadAnalysis 5 database. However, a comparison of the observed limits in both cases reveals that there is no significant deviation when updating to the larger luminosity. We therefore follow the procedure outlined in Ref. [57] to estimate the full CMS run 2 sensitivity from the CMS-SUS-17-001 one. We recalculate CL s exclusions by extrapolating the background and keeping its relative uncertainty constant, and by assuming that the new number of background events to be equal to the number of observed events.
The ATLAS CONF 2019 040 analysis focusing on a multi-jet final state yields relevant bounds from the ttSS channel. Those bounds are shown in figure 8 (left) in orange in the (m T , m S ) mass plane. In such a plane, the grey area is kinematically forbidden as the DM candidate S is required to be lighter than the mediator T . The lighter and darker orange regions correspond to the 68% and 95% confidence level (CL) excluded regions respectively. As there is no information on the correlations between the different signal regions of the ATLAS analysis, we derive our exclusion bounds by solely considering the most constraining of all signal regions of the analysis. For light dark matter, mediator masses ranging up to 1.25 TeV are excluded. Such a strong bound originates from the associated split spectrum configuration that gives rise to a fair number of hard final-state jets, produced in association with a lot of missing transverse energy. Such a topology being the primary target of the ATLAS CONF 2019 040 analysis, we end up in a situation where the sensitivity is maximised. With the dark matter mass increasing, the average transverse momentum of the jet decreases once we enforce the mediator mass to be not too large so that the NLO QCD production rate stays reachable at the LHC run 2 [39]. In addition, the amount of missing transverse energy decreases accordingly, so that the sensitivity drops when the mediator mass is relatively large (too small fiducial cross section when the p T requirements on the jets are imposed) and small (too compressed spectrum) for a given m S value. For instance, for m S ∼ 500 GeV, we observe that mediator masses in the [800, 1100] GeV range are excluded at 95% CL. Furthermore, for m S 550 GeV we lose all sensitivity. Those bounds significantly improve previous collider limits on dark matter models with a top-philic vector-like portal that are associated with the multi-jet plus / E T search channel. The improvement corresponds to a factor of about 1.3 on the mediator mass, and to a factor of about 2 on the dark matter mass after a comparison with the bounds derived on the basis of early LHC run 2 results [9].
We now turn to bounds extracted from searches for the BSM production of top antitop pairs in association with missing transverse energy. Rescaling to the full LHC run 2 luminosity the limits derived from the CMS-SUS-17-001 analysis, we present the resulting bounds as the lighter and darker blue regions of the left panel of figure 8. These areas are excluded at 95% and 68% CL respectively, and we use once again the best signal region of the analysis to conservatively estimate our bounds. This analysis targeting precisely the main collider signature of the considered model (pp → ttSS), we can expect quite a high sensitivity to the signal. We indeed find that dark matter masses as high as about 700 GeV are reached (for a heavy mediator of about 1 TeV, so that the spectrum is not too compressed), which extends the reach of the multi-jet plus / E T search. In addition, the tt plus / E T bounds also complement the multi-jet ones in the more compressed regime, for mediator masses lying in the [400, 1000] GeV mass window (and for dark matter being correspondingly at least 40% lighter). For very compressed spectra, the final-state objects are not hard enough in general, so that the analysis loses sensitivity, similarly to the the multi-jet plus / E T case. Moreover, larger mediator masses are also hard to probe due to the steeply falling NLO QCD production cross section. In order to quantify the improvements of the sensitivity relatively to the early LHC run 2 results, figure 8 includes as a blue dashed line the early run 2 exclusion at 95% CL that relies on a luminosity of 35.9 fb −1 . This shows that a factor of 4 in luminosity allows for stronger bounds when the mediator is heavy and the spectrum is very split, the m T lower limit increasing by about 200 GeV, for almost no change in the compressed regime.
Finally, in the right panel of figure 8, we extrapolate our results to 3 ab −1 , which corresponds to the expected luminosity of the high-luminosity operations of the LHC. This extrapolation follows the strategy outlined above. The parameter space area that is covered extends by about 10%-15%, both in terms of dark matter and mediator masses when the mass spectrum is split. On the contrary, the compressed regime shows almost no improvement, as was already the case for the comparison of the bounds obtained when using a luminosity of 35.9 fb −1 and 139 fb −1 . For these compressed scenarios, as will be discussed in the next section, cosmological probes are however in order to probe the model. A large fraction of the parameter space featuring cosmological properties in agreement with current data will however stay un-probed for the next decades. When both the dark matter and the mediator lie in the TeV or multi-TeV regime, there is indeed currently no sensitivity, either cosmologically or from colliders. Such a region being out of reach of the LHC and any planned dark matter direct or indirect detection experiment, the most fruitful strategy may be to rely on a future proton-proton collider option that would run at 100 TeV. The estimation of the corresponding reach is left for future work, where the range of validity of the effective operator should be treated with care.

Summary and conclusions
In this work we have examined a top-philic scalar dark matter scenario as could emerge from composite Higgs models. Our study relies on an ad hoc simplified model at the TeV scale. This model features a heavy top-philic DM candidate S and a heavier vector-like fermion T mediating the interactions of the dark matter with the top sector. In contrast to previous investigations, we have included not only a Yukawa coupling of the form ST t, but also a dimension-five contact operator ttSS, as both are generally predicted in composite setups. Focusing on scenarios in which the heavy vector-like top partner has a mass comparable to the DM mass, we have investigated the parameter space which yields the correct relic density numerically and semi-analytically (see section 2 and appendix A for details). We have in particular examined the interplay between the contributions at dimension-four (proportional to the ST t Yukawa coupling) and those at dimension-five (proportional to the ttSS Wilson coefficient). In investigating the direct detection constraints (see section 3), we have shown that the inclusion of the dimension-five operator impacts the determination of the viable region of the scalar top-philic DM parameter space in a significant way. DM is hence allowed to be as heavy as 3 TeV in many different configurations, the corresponding direct detection cross sections being below the current XENON 1T bounds and either above (with future detection prospects) or below (and thus hard to probe) the neutrino floor. We have next studied constraints and expected bounds arising indirect DM detection in section 4, emphasising again the important role played by the ttSS operator. Finally, current and projected LHC bounds, that are in contrast agnostic of the considered dimension-five operator, have been examined in section 5.
We summarise our results in figure 9 which shows the various bounds and the still allowed parameter space regions. The exclusions are shown in a (m S , r = m T /m S −1) plane, or in other words in a plane with the dark matter mass m S and the spectrum compression factor r as x and y axes respectively. We consider three setups. In the first one (top panel of the figure), we consider, for each pair of S and T masses, the largest possible value for the Wilson coefficient C/Λ for which there exists aỹ t value leading to the right relic density. In the second and third considered configurations, we choose instead a scenario in which |C/Λ| is minimum, theỹ t values being derived to reproduce the right relic density. We distinguish scenarios featuring a negative C/Λ value (second panel of the figure) and a positive one (third panel of the figure). We can immediately observe the dependence of the cosmological bounds on the C value, the collider bounds solely depending in contrast on the new particle masses. A larger C value implies weaker cosmological bounds, whilst a minimum C value reduces the impact of the direct detection bounds and increases the one of the indirect detection ones in a complementary manner. In addition, direct and indirect detection probes are the only ones relevant to enter the compressed region in which r is small. Colliders have no or very little sensitivity in this regime. On the contrary, future colliders are the only way to access the so-far allowed large-mass region of the parameter space. Without such machines (whose sensitivity estimation lies beyond the scope of this work), scenarios featuring a dark matter mass larger than about 700 GeV may stay unreachable, although they are fully viable in the light of reproducing the DM relic density as observed by Planck.
in Lyon. LM is supported by the UJ GES 4IR initiative, and thanks Campus France for support under the Eiffel programme. TF's work is supported by IBS under the project code IBS-R018-D1.

A Details on the relic density fit
In the following we go into further detail on the semi-analytical fit of the curve which relates the parametersỹ t and C/Λ in producing the correct relic density. This curve depends on both m S and m T for its shape, so that imposing that the relic density matches Planck data leads to where the function f has to be determined.

A.1 Without co-annihilations
We illustrate the behaviour that the f function should reproduce in figure 10 for several benchmark configurations. In this figure, the behaviour of the curve is studied for approximately constant r = m T /m S −1 values across a number of m S benchmarks, and we consider scenarios for which co-annihilations are negligible. The treatment of the co-annihilation is left for the next subsection. First, we can notice that the value of m S determines the value of C/Λ for which the dimension-five operator takes over entirely from the Yukawa coupling. In other words, for each DM mass there exists a maximum Yukawa coupling value so that C/Λ has to be constant to reproduce the observed relic density. This is to be expected, as σv dim5 depends only on m S , and not on m T as in eq. (2.5). This dependence of C/Λ on m S can be seen from the green and red lines in figure 10, where the difference in r (or equivalently on m T ) modifies the slopes but not the value of C/Λ which takes over for small Yukawa values. Additionally, the figure shows that for constant r, a modification in m S changes both the maximal C/Λ value and the shape of the curve. It is evident that the addition of the dimension-five operator ttSS allows for a range of Yukawa couplingsỹ t for the ST t operator, where previously (i.e. without adding the ttSS operator to the Lagrangian) only a singleỹ t was allowed to get the correct relic for a given mass point. In particular, we find that the addition of the dimension-five operator allows for viable scenarios featuring a relatively small Yukawa coupling. In this regime, the relic density is entirely driven by the ttSS Lagrangian term. This is the first notable value of interest, where this C/Λ lies between 0.174 and 0.166 TeV −1 . For smaller values of C/Λ, the contribution from the dimension-five operator to the cross section is lower, and the relic is brought back to Planck's value by the contributions involving the ST t operator. In short, the addition of the dimension-five operator extends the viable part of the parameter space, allowing smaller values of the Yukawa couplingỹ t which were previously forbidden.
For r values such as the one used in figure 10, we obtain the fitting functions shown in eq. (2.11). ] Figure 10. A comparison of theỹ t versus C/Λ curve yielding the correct relic density for a number of mass configurations. Although a consistent shape is observable, the gradient and shift of the curve is clearly dependent on the mass parameters.

A.2 Parametrising the shift due to co-annihilations
While for larger values of r we are able to straightforwardly fit directly C/Λ from eq. (2.11), for smaller values of r (where S and T are closer in mass), we find a deviation from the fit. A small shift along the x-axis (i.e. a shift inỹ t ) is observed between the predictions from our scan and the fitting function. This phenomenon is visible in figure 11, where the red line (the initial fit) deviates from the data points (i.e. our numerical predictions), shown over varying r for an illustrative scenario with m S = 682.6 GeV. For smaller values of r (where the "smallness" that is relevant depends on the masses of the particles), we find that the predictions are subject to a shift inỹ t . We estimate this shift for the considered scenarios by the yellow lines of figure 11, in which the fit in shifted by some constant amount in order to agree with the relic density predictions. The deviation between the predictions and the fitting function can be modelled simply with good agreement, as the shift follows an exponential growth as r gets smaller. The amount by which the function must be shifted hence gets exponentially larger as r grows smaller. This is further illustrated in the lower right panel of figure 11. This feature can be understood by examining the impact of co-annihilation processes on the relic density, which is relevant only when the DM candidate S is close in mass with another resonance (in this case, the mediator T ). In this scenario, the relic abundance is driven not only by self-annihilation, but also by co-annihilations between S and T , which leads to the annihilation cross section no longer being given by the simple eq. (2.9).
The calculation of the relic density must be modified to the co-annihilation case in a generalised fashion [58][59][60]. The Boltzmann equation (2.6) is generalised to a set of coupled equations governing the evolution of the S and T states through the universe's history. Focusing on the dark matter (co-)annihilation cross section only, we have The fitting function f from eq. (2.11) fits well for more separated masses of S and T (bottom left figure), but when the mass difference is small (figures of the upper row) the function deviates from the observed values. The function can be shifted in the positive x direction (that is, inỹ t ) in order to re-establish agreement. In the bottom right plot, the deviation is fitted to a decreasing exponential, where the largest shifts inỹ t are necessary for lower r, and the shift becomes negligible above r = 0.8. All points yield the correct relic density.
where x = m S /T * with T * being the temperature, and where g S and g T stand for the DM and mediator number of internal degrees of freedom, and g ef f for the effective number of internal degrees of freedom. Moreover, σ SS and σ ST respectively correspond to the annihilation (SS → tt) and co-annihilation (ST → t ( * ) → X) cross sections. It is then apparent that we can expect the deviation from the fit to be exponentially larger for smaller r, and we can approximate this deviation using an exponential function.
Across the range of allowed masses for the DM candidate, two separate regimes are observed in the behaviour of the interplay of the different contributions to the annihilation cross section. Below m S ≈ 1 TeV, the NLO effects are not yet dominant, and σ N LO /σ vqq ≈ 1. For higher mass regions, the NLO cross section is influenced by the VIB contributions, as apparent in the example benchmarks in figure 2. This motivates distinguishing two mass regimes, above and below about 1 TeV. At each m S benchmark we hence find a smooth pattern of shifts across r which can be fitted to an exponential function. Repeating the fit done in figure 11 across a number of benchmarks, the function which determines the shift where k is an unconstrained parameter of dimension 1 allowed to vary across the benchmarks (but constant for a given m S ). The constants are determined by the fit. In figure 12 (left), we show these exponential functions which map the shifts for a number of m S benchmarks. Finally, we would like to be able to estimate the value of the parameter k in the exponential shift function, as it is specific to the mass m S . We find that the values of this also may be fitted to an exponential function fully determined by m S , as displayed in figure 12. We obtain k(m S ) = (1.9 × 10 −4 ) 6.2 × 10 8 m S /Λ m S ≤ 1.2 TeV for Λ = 3.5 TeV. The value of the dimension-five coupling contributing to the relic density for a given benchmark (m S , m T ) can then be fully determined by extending eq. (2.11) to include the shift, as previously quoted in eq. (2.13).