Top-philic dark matter within and beyond the WIMP paradigm

We present a comprehensive analysis of top-philic Majorana dark matter that interacts via a colored t-channel mediator. Despite the simplicity of the model -- introducing three parameters only -- it provides an extremely rich phenomenology allowing us to accommodate the relic density for a large range of coupling strengths spanning over six orders of magnitude. This model features all `exceptional' mechanisms for dark matter freeze-out, including the recently discovered conversion-driven freeze-out mode, with interesting signatures of long-lived colored particles at colliders. We constrain the cosmologically allowed parameter space with current experimental limits from direct, indirect and collider searches, with special emphasis on light dark matter below the top mass. In particular, we explore the interplay between limits from Xenon1T, Fermi-LAT and AMS-02 as well as limits from stop, monojet and Higgs invisible decay searches at the LHC. We find that several blind spots for light dark matter evade current constraints. The region in parameter space where the relic density is set by the mechanism of conversion-driven freeze-out can be conclusively tested by R-hadron searches at the LHC with 300\,fb$^{-1}$.


I. INTRODUCTION
Astrophysical and cosmological probes continue to consolidate our knowledge of the gravitational impact of dark matter (DM) (see e.g. [1,2]). In order to identify its nature and pinpoint its interactions with the standard model (SM) it is required to explore the cosmologically viable parameter space of models incorporating a DM candidate and confront it with experimental constraints exploiting the large amount of results from DM searches (see e.g. [3,4] for reviews). A common framework to perform such a study in a bottom-up approach are simplified models (see e.g. [5,6] for reviews and references therein) which are assumed to describe the physics at the phenomenologically relevant scales of a (possibly more complicated) UV-complete theory to good approximation.
In theories of new physics the top quark and its couplings play a special role due to a possible link to Higgs physics: Models beyond the SM alleviating the gauge hierarchy problem frequently invoke top partners, such as the stop within the minimal supersymmetric standard model (MSSM). Here we study a DM simplified model with a Majorana fermion DM candidate and top-philic t-channel mediator. 1 Despite the simplicity of the model it provides an extremely rich phenomenology. For coupling strengths in the ballpark of the SM gauge couplings the relic density can be generated by DM freeze-out with or without strong coannihilation effects depending on the mass split-ting between the DM and the mediator. In this region DM shares the properties of a weakly interacting massive particle (WIMP). This region has been studied in the context of supersymmetry [22][23][24][25][26][27] or more generally considering a free coupling strength [28,29]. Here, we confront the model with updated constraints from the LHC as well as direct and indirect detection. We include the region of DM masses below the top mass where loopinduced and three-body final state annihilation channels are important. Furthermore, we consider direct detection signals as well as DM production at the LHC via one-loop processes which can be probed by invisible Higgs decay and monojet searches.
For much smaller couplings DM freeze-out has to be revised since the commonly made assumption of chemical equilibrium between DM and the mediator cannot be maintained during freeze-out. This leads to a phenomenologically distinct variant of DM genesis where the relic density is primarily determined by the rate of conversion processes between DM and the mediator [30,31]. The phenomenological consequences are striking: The parameter region accommodating the measured relic density via conversion-driven freeze-out cannot be probed by conventional WIMP searches but predicts new signatures of long-lived particles at colliders [30]. Despite the small couplings conversion-driven freeze-out allows for thermalization of DM and hence washes out any dependence on the thermal history prior to freeze-out -an appealing feature that is maintained from the WIMP parameter region. With respect to the analysis in Ref. [30] which considers the same simplified model but for a bottompartner mediator several quantitative differences arise. Due to the non-negligible mass of the top, the decay rate of the mediator is kinematically suppressed and conversions are dominated solely by (2 → 2 and 2 → 3) scatterings. As the decay becomes efficient only well after freeze-out the scenario is subject to constraints from Big Bang nucleosynthesis (BBN). Furthermore, the mediator becomes stable on typical timescales for traversing an LHC detector.
In this study we present a comprehensive analysis considering both regions. The structure of the paper is as follows: After introducing the model considered in this work in Sec. II we discuss the cosmologically viable parameter space in Sec. III. We categorize the parameter space in the WIMP and conversion-driven freeze-out scenario presented in Secs. III A and III B, respectively. The latter section includes a detailed discussion of the underlying Boltzmann equations in the out-of-chemicalequilibrium regime. In Sec. IV we confront the cosmologically preferred parameter space with experimental constraints from direct detection, indirect detection, various collider searches and BBN constraints. We conclude in Sec. V.

II. SIMPLIFIED TOP-PHILIC MODEL
We consider a simplified model containing a neutral Majorana fermion χ that transforms as a singlet under the SM gauge groups and a colored scalar particlet with gauge quantum numbers identical to the right-handed top quark. We furthermore assume a Z 2 symmetry under which all SM particles are even while χ → −χ andt → −t are odd. Under these assumptions the Majorana fermion χ is absolutely stable for m χ < mt and provides a viable DM candidate. The interactions of these particles with the SM are described by the Lagrangian where D µ is the usual covariant derivative and t the top quark Dirac field. The coupling λ χ characterizes the coupling strength of the DM particle with the SM, being mediated by the colored scalart. The simplified model is characterized by three parameters, the masses m χ , mt and the coupling λ χ . The particle content and interaction terms can be viewed as being part of the stop-neutralino sector of the MSSM. More specifically,t corresponds to the righthanded stop and χ to the bino in the supersymmetric context. The coupling constant is then fixed to 3 cos θ W ≈ 0.33. However, the simplified model can also be considered as the low-energy limit of nonsupersymmetric extensions of the SM. Alternatively, if supersymmetry is realized in nature, it could be nonminimal, i.e. described by a particle content beyond the MSSM. For example, χ = sin(θ)B 0 + cos(θ)S can be a mixture of the binoB 0 and the fermionic componentS of an additional supersymmetric multiplet that is a SM singlet. In this case the coupling λ χ = sin(θ)λ MSSM χ will be reduced compared to the MSSM value by the mixing angle [32]. In the following, we will be ignorant about the embedding of the simplified model within extensions of the SM and treat the coupling λ χ as a free parameter.
Note that the gauge and Z 2 symmetries allow an additional renormalizable interactiont †t H † H to the Higgs field. It would lead to a correction of the t t * annihilation cross section, contributing to the coannihilation rate, of the scattering rate χN → χN off nuclei via Higgs exchange, relevant for direct detection, and of the loopinduced annihilation rate via a Higgs in the s-channel, that can become important for m χ ∼ m h /2. If the corresponding coupling is well below unity, its effect is subleading compared to the processes mediated by strong and top-Yukawa interactions, respectively [28]. In the following we assume that this is the case.
In addition, a flavor-violating coupling oft to righthanded charm or up-type quarks can be considered. If present, it could potentially have a sizeable effect on the lifetime oft for small mass splittings ∆m = mt − m χ m t . Even if we impose that flavor-violating couplings vanish at a certain energy scale µ 0 , they are generated by renormalization group (RG) running at a different scale µ. For example, the RG-induced coupling of the same form as in (1) to charm quarks is given by [33] λ c where the numerical estimate corresponds to RG running between the scale of grand unification and the electroweak scale. We find that the impact of this RGinduced flavor-violating coupling on the decay of the mediator can safely be neglected for ∆m > ∼ 10 GeV. For even smaller mass splitting a certain degree of tuning would be required to suppress flavor-violating decays. For the purpose of this work we assume in the following that the relevant interactions are captured by the Lagrangian given in Eq. (1).

III. DARK MATTER FREEZE-OUT
The simplified model encompasses two regions in parameter space for which the processes that are responsible for setting the DM abundance are qualitatively different. First, there is a portion of parameter space where either DM annihilation or coannihilation processes involving the mediator t govern the relic density, and to which we refer as "WIMP region". Second, for small enough value of the mass splitting ∆m = m t − m χ and m χ < ∼ 2 TeV, the dark matter density is set by conversiondriven freeze-out. In this region of parameter space, the commonly adopted assumption that conversions between t and χ are in equilibrium during dark matter freeze-out breaks down. We discuss both regions in the following two subsections, respectively.

A. WIMP region
In the WIMP region the relic density is set by the usual freeze-out mechanism. Depending on the relative size of m t and m χ , also coannihilation processes have to be taken into account. Within the WIMP region the coupling λ χ is large enough in order to guarantee chemical equilibrium between DM and the coannihilation partner such that the relic density is determined by the effective, thermally averaged annihilation cross section following the common coannihilation scenario [34], where n eq = i n eq i and n eq i = T /(2π 2 ) g i m i K 2 (m i /T ) with g χ = 2 and g t = g * t = 3. Let us now discuss the dominant annihilation channels within various parts of the WIMP region. For m χ > m t , the annihilation channel χχ → tt is kinematically allowed. In addition, the coannihilation channels χ t → tg and t t * → gg give a sizeable contribution to σv eff if m t /m χ < ∼ 1.3. The corresponding cross sections scale as λ 4 χ , λ 2 χ g 2 s , and g 4 s , respectively, where g s is the strong coupling constant evaluated at a scale µ ∼ m χ . Note that the t t * annihilation cross section is independent of the value of λ χ , such that it dominates for small values of λ χ as long as the assumption of chemical equilibrium is justified. We will discuss the breakdown of this assumption in the next subsection.
For m χ < m t annihilation into a pair of top quarks is kinematically disfavored, leading to a strong reduction of the DM annihilation cross section. For the computation of the relic density we include the loop-induced annihilation χχ → gg into a pair of gluons [35] as well as the 2 → 3 process χχ → W bt [36][37][38] (here and below W bt stands for the sum of W + bt and W −b t). In addition, we include the channels χχ → h ( * ) → bb, W W * , . . . arising from the loop-induced Higgs-χ-χ coupling [39], which is resonantly enhanced for m χ ∼ m h /2, where m h = 125 GeV (see below for details).
For m χ + m t < m t also the coannihilation process χ t → tg becomes kinematically forbidden. In this case the annihilation channel χ t → W b, involving the weak instead of the strong coupling, becomes important. The top-quark in the s-channel leads to a resonance for m χ + m t m t , such that this process has a large impact for very small masses.
The relic density within the WIMP region is computed based on a modified version of micrOMEGAs 4.3.1 [40]. Apart from the 2 → 2 processes that are included by default, we added the loop-induced process χχ → gg, χχ → h ( * ) → bb, W W * , . . . and the 2 → 3 channel χχ → W bt. In addition we also include the effect of Sommerfeld enhancement as described in Appendix B of [30]. Possible further refinements include bound-state effects [24,[41][42][43][44][45] as well as next-to-leading order (NLO) corrections to coannihilation rates [46,47] within QCD, which are, however, beyond the scope of this work. The full annihilation cross section for χχ → gg has been extracted from Ref. [35] and is too lengthy to report here. In the limit m χ m t , m t it is given by where f (r) ≡ 1 2(r − 1) 4 r 3 − 6r 2 + 3r + 6r ln(r) + 2 . (5) For m t m χ , m t one finds where δ = m t /m χ − 1.
The cross section for 2 → 3 annihilation is given by where µ i = m 2 i /m 2 χ , y = E t /m χ and γ t = Γ t→W b /m χ are the energy of the final-state top and the top width normalized to the DM mass, respectively, and λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2zx is the Källén function. The integration boundaries are given by y min = m t /m χ and The process is taken into account below the tt threshold.
For each pair of masses (m χ , m t ) we fix the coupling λ χ such that the relic density resulting from freeze-out matches the measured DM density Ωh 2 = 0.1198 ± 0.0015 [1]. In Fig. 1 we show the resulting contour lines of constant coupling λ χ , where we use the DM mass and the mass splitting ∆m = m t − m χ as independent parameters. We also indicate explicitly the contour for which λ χ matches the bino-stop-top coupling within the MSSM. If we restrict the coupling to be less than 4π, the relic density exceeds the measured value within the grey-shaded region, and we therefore disregard it in the following. Below the thick black line coannihilations are so efficient that the relic density resulting from the freeze-out computation described above would lie below the measured value. This parameter domain is discussed in detail below. The remaining part of parameter space corresponds to the "WIMP region".
The kinematic threshold for tt annihilation is clearly visible in the contours shown in Fig. 1, and leads to the sharp drop for m χ ∼ m t . For m χ < ∼ m t the annihilation channel χχ → W bt as well as the Boltzmann tail of the DM distribution allowing for χχ → tt yield a sizeable contribution, and slightly smoothen the step-like behavior of the contour lines. Coannihilations start to play a role for ∆m < ∼ m χ /10, and allow for larger DM masses for a given coupling. Additionally, for very small masses, the contours feature a 'spike' at m χ ∼ m h /2 as well as a 'bump' for m χ + m t m t . The 'spike' can be explained by the Higgs resonance in the loop-induced annihilation channels χχ → h ( * ) → bb, W W * , . . . , and the 'bump' at slightly higher mass is related to a top resonance in the coannihilation process χ t → W b.

B. Conversion-driven freeze-out solutions
As mentioned above, up to a DM mass of around 2 TeV we encounter a region in parameter space with small ∆m where the effective, thermally averaged cross section for mediator-pair annihilation alone -which is fixed for a given DM mass and mass splitting -is so large that one undershoots the measured relic density, seemingly regardless of the coupling λ χ . However, this conclusion hinges on the assumption of chemical equilibrium between DM and the mediator, i.e. the condition n χ /n eq χ = n t /n eq t , which does not hold once λ χ decreases beyond a certain value. In fact, dropping this assumption one can find solutions with small λ χ where the relic density is governed by the mechanism of conversion-driven freeze-out [30]. In the following we will first outline the computational steps of the relic density calculation before we discuss the phenomenological aspects in the corresponding region in parameter space.

Boltzmann equation and conversion rate
In the absence of chemical equilibrium between DM and the mediator the computation of the relic density requires us to solve the coupled set of Boltzmann equations for the respective abundances [30,51,52], where Y = n/s is the comoving number density, s the entropy density and x = m χ /T is the temperature parameter. Y t represents the summed abundance of the mediator and its anti-particle. This leads to the factor 1/2 in Eq. (12), since the cross sections are averaged over initial state degrees of freedom.
In addition to the terms accounting for annihilation and coannihilation displayed in the first lines of Eqs. (11) and (12) three further terms occur in the second lines in both equations that account for conversion processes. The first of these terms includes conversion via scattering processes. As will be discussed below, both 2 → 2 as well as 2 → 3 and 3 → 2 processes need to be included due to the Boltzmann suppression of W and t in 2 → 2 scatterings at low temperatures (see Fig. 2 for an illustrative example). The respective rate is given by where and k, l, m denote SM particles (see below). Γ χ→ t is understood to contain the conversion into both the mediator and its anti-particle which leads to the factor of two in front of the sums in Eqs. (14) and (15). In the last line we used Γ χ→ t Y eq χ = Γ t→χ Y eq t to rewrite the conversion rate such that it contains two-body initial states only, which turns out to be more convenient for numerical evaluation. In the following we refer to the sum of the second and third line as 2 → 3 conversion processes for brevity. Neglecting quantum statistical factors and assuming thermal momentum distributions the thermally 2. Examples for diagrams contributing to 2 ↔ 3 conversion processes which are taken into account below the kinematic threshold of the 2 → 2 scattering on the left.
averaged cross sections are given by [34] σ ij v n eq i n eq j = where g i are the internal degrees of freedom of species i, p ij denotes the absolute value of the three momentum of the initial state particles i, j in the center-of-mass frame and K n denotes a modified Bessel function of the second kind. The middle terms in the second lines of Eqs. (11) and (12) account for the conversion induced by decay and inverse decay. The thermally averaged decay rate reads where γ is the Lorentz factor. For small mass splitting ∆m the leading contribution is the 4-body decay t → χbf f , where f, f are light SM fermions. Finally, the last terms in the second lines of Eqs. (11) and (12) take into account the scattering processes within the odd-sector. We now describe the various conversion processes that are taken into account in detail. We compute the squared matrix elements |M | 2 for all 2 → 2 processes with CalcHEP [53] including all diagrams that are allowed at tree-level. For the computation of the 4-body decay t → χbf f we use Madgraph5_aMC@NLO [54]. Furthermore, we include the following 2 → 3 scattering pro-cesses: where on the right-hand side we show the 2 → 3 processes that supersede the 2 → 2 shown on the left-hand side below the respective kinematic threshold. Figure 2 shows examples for the respective diagrams. The fermions f, f are all possible final states the W -boson is allowed to decay into (i.e. light quarks, neutrinos and charged leptons). We calculate the scattering cross sections as function of s with Madgraph5_aMC@NLO in the direction 2 → 3 as discussed above. Since the initial state particles in the various processes have different masses we integrate each contribution separately up to the threshold of the corresponding 2 → 2 processes using Eq. (17) before combining the rates. We do not include the contributions (similar processes exist for γ, Z, H instead of g) as the rates for these processes are of the order of O(α s ) corrections to the original 2 → 2 processes from Eq. (19). Their consistent inclusion would require a full NLO computation of the conversion rate, including also virtual corrections, which is beyond the scope of this work. For a detailed discussion of possible deviations from kinetic equilibrium we refer to Appendix C of Ref. [30], which were found to have only a minor impact on the final DM abundance.

Phenomenology of conversion-driven freeze-out
In the following we discuss the properties of the solutions of Eqs. (11) and (12) taking the deviation from chemical equilibrium in the odd-sector into account. We focus in particular on the differences compared to the bottom-partner mediator considered in [30].
In the left panel of Fig. 3 we show the time evolution of the DM and mediator abundance for an exemplary parameter point m χ = 1400 GeV, m t = 1420 GeV, λ χ = 4.6 × 10 −6 . Here the coupling λ χ has been adjusted such that the final DM abundance matches the observed value. The corresponding (co)annihilation and conversion rates are shown in the right panel of Fig. 3. Shortly after the dominating annihilation rate of the mediator drops out of equilibrium (blue solid line), also the conversions freeze out (red dashed line). This leads to a reduced conversion χ → t such that the DM abundance is not depleted as strongly as it would be when conversions are in equilibrium.
The large top mass leads to several qualitative (and also quantitative) differences compared to the coupling to bottom quarks considered in Ref. [30]. In particular m t is larger than ∆m in the region that potentially allows for conversion-driven freeze-out. Accordingly, the 2-body decay and, in fact, also the three-body decay of the mediator is kinematically forbidden rendering the decay rates to be small compared to the scattering rates (solid and dashed red lines in the right panel of Fig. 3). For this reason the decay rate becomes efficient only well after freeze-out so that freeze-out and decay take place separated from each other. This can be seen in the time evolution of the abundances shown in the left panel of Fig. 3.
Another difference to the case of a bottom partner mediator is the necessity to include 2 → 3 conversion processes. In particular in the region of small DM masses we encounter a larger difference between the relic density prediction with or without 2 → 3 conversion processes. The difference can be seen in Fig. 4 for m χ = 100 GeV, m t = 110 GeV where the solution for λ χ that accommodates the measured relic density differs by more than an order of magnitude. This is due to the strongly suppressed top and W ± abundances at temperatures relevant for freeze-out for DM masses around or below the top mass. For larger DM masses the importance of 2 → 3 conversion processes sets in at later times where the respective rates are, however, already not efficient anymore, cf. the red dashed curves at x > ∼ 200 in the right panel of Fig. 3. Hence, for large m χ the relic density is governed by the 2 → 2 processes, which are significantly larger. This explains why the coupling λ χ required to explain the observed DM abundance depends strongly on the DM mass and ranges from ∼ 10 −3 around m χ = 100 GeV to ∼ 10 −6 for DM masses around 1.5 TeV. The contours of λ χ are shown in Figs. 1 and 9 below the thick black line, which corresponds to the region in parameter space where conversion-driven freeze-out is important.
A peculiarity of the present model is the appearance of a small region in parameter space providing solutions of Ωh 2 = 0.12 for three different values of the coupling λ χ , for given DM and mediator mass. The region where this occurs lies within a thin band just outside the boundary of the conversion-driven freeze-out region (as indicated by the thick black curve in Fig. 1) towards large m χ . The corresponding functional dependence of Ωh 2 on λ χ is shown in the left panel of Fig. 5. While the blue curve (corresponding to DM and mediator masses just inside the boundary of the conversion-driven freezeout region) allows only for one solution, the red curve (a point somewhat outside this boundary) exhibits three solution due to a local minimum of the function around λ χ 10 −5 . This minimum appears somewhat counterintuitive at first sight. To understand its origin we consider the ratio of the deviations from thermal equilibrium of χ and t, i.e. the quantity Y χ /Y eq The right panel of Fig. 5 shows the evolution of this quantity with x for the three parameter points indicated in the left panel. For typical WIMP freeze-out (point C) chemical equilibrium is maintained and this quantity is equal to unity. In contrast, for point A, which lies within the region of conversion-driven freeze-out, this quantity deviates significantly from one. However, the deviation occurs in both directions, depending on x: At early times up to x 100 DM is relatively over-abundance with respect to the mediator. This stems from the fact that the mediator annihilates efficiently while the conversion ratesthat are responsible for reducing the DM abundanceare on the edge of being efficient. That is, the DM abun-  dance is not able to entirely follow the fast reduction of the mediator abundance. At late times, however, the mediator is relatively over-abundant with respect to the value in chemical equilibrium. This occurs at x > ∼ 100 where T ∼ ∆m. At around this temperature conversions start to prefer the direction t → χ. If these conversions are less efficient than needed to maintain chemical equilibrium it leads to the observed relative over-abundance of the mediators. Finally, for point B, which is close to the observed minimum, the transition to a relative overabundance of the mediators appear earlier while late annihilations are still not entirely inefficient. The larger mediator abundance leads to a slightly larger annihilation rate which reduces the overall abundance in the dark sector. For a parameter point with a somewhat larger coupling (towards the plateau) the conversion rates are larger reducing the relative over-abundance of the mediator and hence reducing the efficient mediator annihilation. The appearance of the three solutions can also be seen in the left panel of Fig. 6, by the slight bending of the vertical part of the relic density line for m χ 1.5TeV.
Due to the small coupling λ χ one may wonder whether the final DM abundance depends on the initial condition at small x. We checked that this is not the case by varying the initial DM abundance over two orders of magnitude with respect to the equilibrium abundance, which has only a negligible impact on the final value. The reason for this is that conversion processes are more efficient at small x (see right panel of Fig. 3), thereby erasing any dependence on the initial condition. This feature distinguishes conversion-driven freeze-out from other beyond-WIMP scenarios, for which dark matter is never in thermal equilibrium, like the superWIMP scenario [55]. Note that this scenario could also be realized in the present model, when considering an even smaller coupling λ χ . In this case the mediator freezes out decoupled from dark matter and decays well after its freeze-out. Assuming a zero initial DM abundance the corresponding allowed parameter space (Ωh 2 = 0.12) is constraint to the black dotted curve shown in Fig. 1 for illustration. 2 However, we do not further consider this case here.

IV. EXPERIMENTAL CONSTRAINTS
In this section we confront the cosmologically allowed parameter space with a wide range of experimental constraints. We consider direct and indirect detection, collider searches for missing energy and long-lived particles as well as constraints from BBN.

A. Direct detection
In the considered model spin-independent DM nucleon scattering is induced by two processes. First, via an effective DM-Higgs coupling induced by triangle diagrams with top quarks and mediators in the loop [39]. Second, a coupling to the gluon content in the nucleus is induced through box diagrams again with top quarks and mediators in the loop [56]. Loop-induced couplings to the Z-boson are suppressed due to the Majorana nature of the DM particle. 2 We plot the curve for which (Ωh 2 )χ = mχ/m t (Ωh 2 ) t = 0.12 where (Ωh 2 ) t is the freeze-out abundance of the mediator in the absence of any coupling to DM. For considerable mass splittings where the 2-body decay of the mediator is open, consistency with BBN (cf. Sec. IV F) can easily be achieved (e.g. τ t < 1 s requires λχ > ∼ 10 −12 for which Γconv/H 1 until well after mediator freeze-out).
The spin-independent cross section for elastic χnucleon collisions is given by [57] where m N is the nucleon mass, and (22) The first term arises from Higgs exchange, with m h = 125 GeV and v = 246 GeV, and involves the hχχ coupling g hχχ given in Eq. (8), which is generated by top/ t loops [28,39]. For the direct detection cross section one needs to evaluate g hχχ in the limit of low momentum transfer s → 0 (see Sec. III A). In this limit Higgs exchange gives rise to an effective dimension-sixχχqq interaction. For the nuclear parameters f N T q (which include the quark masses encoding the coupling to the Higgs boson) we use the values reported in [58] (see also [59]), and f N T G = 1 − q=u,d,s f N T q . The second and third terms in Eq. (22) arise from loopinduced couplings to gluons [56,60] where we use the loop functions I j (m t , m t , m χ ) taken from [60]. The first contribution corresponds to an effective dimension-sevenχχG aµν G aµν interaction, where G aµν is the gluon field strength tensor, and g G is related to the coefficient of the gluon twist-2 operator, which contains additional derivatives acting on χ [60]. Note that b and g G are well-behaved for m t → m t + m χ . For the gluon twist-2 nuclear parameter we use G 2 = 0.48 [61]. For small DM and mediator masses the contribution of the loop-induced coupling b to gluons dominates over the Higgs-exchange contribution. For illustration we provide analytic expressions in the heavy-top limit, where y 2 t = 2m 2 t /v 2 . When increasing the mediator mass the Higgs-exchange contribution becomes more important and ultimately dominates. Due to the relative sign difference both contributions can interfere destructively and lead to a blind spot for certain parameter values. In the heavy-mediator limit one finds, in agreement with [28], The couplings to gluons are suppressed compared to the Higgs-exchange contribution for large mediator masses, as expected. The above expressions for g hχχ agree with Eq. (9) for s → 0 and r = m 2 t /m 2 t → 0 or ∞, respectively. In our numerical analysis of direct detection rates we use the full expressions for the loop-induced couplings for s → 0.
In order to derive 90% C.L. constraints on the model parameter space we compare the cross section to the respective limits from Xenon1T [62]. Furthermore, we show projections for LZ [63]. The results are included in Figs. 6 and 8.
For ∆m = 20 GeV (left panel in Fig. 6), current limits exclude couplings down to λ χ > ∼ 0.7, excluding the thermal relic scenario (green curve) for DM masses below 64 GeV. However, this region is already ruled out by LEP searches (grey shaded area). The destructive interference of the Higgs-exchange and gluon contributions to σ SI also leads to a blind spot for direct detection for m χ 220 GeV, in line with the analytical results discussed above.
For larger values of ∆m the relative importance of the Higgs-exchange contribution becomes larger, and the blind spot correspondingly shifts to smaller DM masses.
For ∆m = m t (right panel in Fig. 6) the cancellation occurs for m χ 40 GeV explaining the decreasing sensitivity towards small masses. The non-observation of a scattering signal by Xenon1T excludes only rather large couplings λ χ > ∼ 3 in this case. Nevertheless, Xenon1T probes the thermal relic scenario (green curve) for DM masses below 135 GeV, with an exception close to the Higgs resonance for m χ < ∼ m h /2. This region is also tested by LHC stop searches (cyan shaded area) as well as invisible Higgs decay for lower masses (light blue shaded area), except for a small range around m χ ∼ 55 GeV (see Sec. IV C for more details).
Projecting onto the surface of the 3-dimensional parameter space that provides a thermal relic (as done in Fig. 8) Xenon1T can exclude a large fraction of the parameter space up to m χ 150 GeV. Despite the raising coupling for large mass splittings Xenon1T slightly loses sensitivity as the scattering cross section is suppressed by the heavy mediator in the loop. In addition, the destructive interference between Higgs exchange and gluonic contributions for particular masses leads to a decrease in sensitivity, which for small mass splitting are approximately given by m χ 225 GeV − 1.1∆m. Finally, direct detection is less sensitive to a thermal relic close to the Higgs resonance at m χ < ∼ m h /2 due to the reduced value of the coupling. The data currently collected by Xenon1T will lead to significant improvements in the near future. Furthermore, the planned experiment LZ is expected to strengthen the sensitivity by a factor of about 2-3 for λ χ , or equivalently 1-2 orders of magnitude for the cross section.

B. Indirect detection
Indirect detection of DM is an important search strategy testing its self-annihilating nature. However, large parts of the parameter space exhibit strong coannihilation effects rendering the direct DM annihilation cross section χχ → SM SM to be relatively small for a thermal relic. Nevertheless, we confront the parameter space with limits from Fermi-LAT γ-ray observations of dwarf spheroidal galaxies, γ-line observations of the Galactic center and cosmic-ray antiproton measurements by AMS-02 as well as projected limits for CTA.
The cross section of DM annihilation today is dominated by four possible channels: (i) The loop-induced χχ → gg channel dominates below and around the W bt threshold, m χ < ∼ (M W + m b + m t )/2, except for a very narrow region where (ii) the (loop-induced) Higgs mediated channel χχ → h → bb, cc, ττ , W W * , ZZ * dominates for m χ m h /2. (iii) The channel χχ → W bt is dominant just below the tt threshold and (iv) χχ → tt dominates above its threshold. Note that for DM masses in the multi-TeV region the annihilation into tt starts to become less important again due to helicity suppression. Therefore annihilation to gg is relevant also for very large DM masses. In addition, 2 → 3 processes (gtt, W bt, Ztt, γtt) can become important in that regime as well [35,38,[64][65][66][67]. Since our main focus is on lower masses we do not consider them here.
We first discuss limits from dwarf spheroidal galaxies. For the prediction of the continuous γ-ray flux we take into account all four channels and sum up their contributions according to their relative weight. For tt, bb, cc, ττ and gg we use the spectrum predictions from [68] which include electroweak corrections from soft and collinear final state radiation. For the three-body final state channel W bt we calculate the spectra using Madgraph5_aMC@NLO [54] and Pythia 8.215 [69]. The spectra for the three-body final state channel W W * and ZZ * are taken from [70]. For the individual cross section prediction we adopt the results discussed in Sec. III A. The predicted energy flux in an energy bin between E min and E max is given by where dN γ /dE γ is the differential photon spectrum per annihilation, J is the J-factor of the considered dwarf and the sum runs over all contributing channels i. We confront the predicted spectra with the Fermi-LAT data [71] using the published likelihoods provided for the individual dwarfs as a function of the energy flux in the considered 24 energy bins. We consider the nine dwarfs with the largest confirmed J-factors as given in [72] and obtain the total likelihood by summing over the individual log-likelihood contributions of all bins for all dwarfs while marginalizing over the J-factor for each dwarf according to its uncertainty provided in [72]. The resulting 95% C.L. limits are shown in Fig. 6 (red solid lines). They reach the thermal relic scenario (green curve with error band) only in a very narrow region around m h /2. In this spot the Higgs mediated annihilation cross section becomes resonantly enhanced for the small DM velocities present today while being less enhanced in the early Universe where the thermal velocity distribution peaks at much larger values. As the width of the resonance is smaller than the widths of the plotted curves in Fig. 6, the limit reduces to a vertical line. It constrains the thermal scenario at m χ m h /2 for mass splittings above ∆m 24 GeV, as indicated by the red arrow in Fig. 8. Below this mass efficient co-annihilation effects allow for a reduction of λ χ and thereby of the indirect detection signal while still allowing us to accommodate the measured relic density.
Next we consider constraints from the cosmic-ray antiproton flux measured by AMS-02 [73] using the results of [74] which provides 95% C.L. limits on σv for various annihilation channels into SM final states. Here we only adopt the limit on tt constraining DM masses above 200 GeV. For smaller masses the cosmic-ray limits are considerably weaker as the analysis exhibits a preference for a DM signal in this region [75]. The results are shown in Fig. 6 (dark red solid lines). The limits come very close to the thermal relic scenario for DM masses between 300 and 400 GeV for ∆m above 30 GeV, i.e. outside the conversion-driven freeze-out region. However, up to relatively large mass splittings this region is also excluded by LHC stop searches. Still, for DM masses above the LHC limit from stop searches (m χ 400 GeV for ∆m = m t ) the antiproton limit places the strongest constraint on the model. Searches for monochromatic γ-lines are a complementary way to probe DM annihilation. In our model annihilation into two monochromatic photons proceeds via the same loop-diagrams as annihilation into gluons. Hence, their cross sections are proportional to each other [76], where we evaluated α s at µ = 300 GeV for the numerical value given. As a consequence γ-line searches can only compete with searches for continuous γ rays for small DM masses for which annihilation into gg dominates. We found that the respective 95% C.L. limits from Fermi-LAT observations of the Galactic center [77] are competitive to the limits from dwarf spheroidal galaxies only for the most aggressive choice of the DM density profile in our Galaxy, namely a generalized Navarro-Frenk-White profile [78] with an inner slope of γ = 1.3. The density profile is, however, subject to large uncertainties. Therefore, and to reduce clutter, we do not show the limit in Fig. 6. Finally, we comment on future prospects for CTA [79] to probe the model. In Fig. 6 (right) we superimpose the optimistic estimate of the projected limit from the observation of the Galactic center presented in [80]. It does not provide sensitivity to the considered model.

C. Stop searches at colliders
At the LHC a large number of searches for neutralinostop simplified models have been performed. As the stop production channel only involves its gauge interactions these searches do not explicitly make reference to the strength of the neutralino-stop coupling. Hence, the results for the decay channels t → t ( * ) χ that do not involve further supersymmetric particles can directly be applied to the model under consideration. 3 However, certain assumptions have to be fulfilled to provide applicability. On the one hand, the width of the (stop-like) mediator has to be sufficiently narrow. Here we require Γ t ≤ 0.2m t as a benchmark value. On the other hand, the mediator decay has to proceed sufficiently promptly in order to match the respective object reconstruction criteria (see below for details).
We consider various 13 TeV analyses, in particular the CMS fully hadronic [81,82], CMS single lepton [83], ATLAS fully hadronic [84], ATLAS single lepton [85] and ATLAS two leptons [86] analyses. For large ∆m the CMS single lepton search [83] provides the strongest bound on the mediator mass reaching m t 1.1 TeV. For smaller ∆m each of the above analyses exhibit certain domains for where it exclusively provides sensitivity. In Fig. 8 we show the 95% C.L. exclusion region from the unification of all analyses listed above (cyan shaded region labeled by 'LHC stop I'). For ∆m < M W + m b only the leptonic searches are relevant, which, however, typically require a small impact parameter of the primary vertex for lepton reconstruction. In order to take this into account we cut the respective limits at a ∆m that corresponds to a proper decay length of 100 µm, see gray short dashed curve that partly marks the lower boundary of the LHC exclusion region in Fig. 8. The gray short dashed curve that partly cuts the exclusion region from above denotes Γ t = 0.2m t .
In addition to the above analysis monojet searches exist that only rely on the missing transverse momentum carried away by χ recoiling against initial state radiation. We consider the 13 TeV ATLAS monojet analysis presented in [87] where it is interpreted within the neutralino-stop simplified model. However, the respective limits are only presented for m t ≥ 250 GeV. We do not assume that the limit extends to smaller masses as the region m t ∼ m t may involve further complications due to similarities of stop and top signals. 4 Hence, the search only constraints a very small fraction of the allowed WIMP parameter space for ∆m < 50 GeV and 210 GeV < m χ < 240 GeV, see Fig. 8 (light gray shaded region denote by 'LHC stop II').
In addition to the LHC limits discussed above limits from LEP provide robust constraints in the region of very small m t . We superimpose the 95% C.L. constraints from data collected by the ALEPH detector presented in [91] which covers the whole range of relevant lifetimes (see [92] for further details). The excluded region is labeled by 'LEP' in Fig. 8.
The exclusion regions are also superimposed in Fig. 6 for ∆m = 20 GeV (left panel) and ∆m = m t (right panel). Note the existence of certain gaps in the low mass region. For ∆m = 20 GeV the stop-searches are not expected to apply to the thermal scenario as a consequence of the required maximum decay length. Although the exact boundary of the region of applicability is subject to some uncertainty, we observe a significant gap between the regions probed by prompt and R-hadron searches (see Sec. IV E for details). For ∆m = m t there exist a small gap for m χ < ∼ 56 GeV.

D. Loop-induced dark matter production
In addition to the searches for mediator production considered in the last section, direct DM production in association with initial state radiation constitutes another search channel. Here we interpret searches for monojet signatures and Higgs invisible decays within the model considering DM masses above and below the Higgs threshold m h /2, respectively.
As the top-content of the proton is negligible, the respective process pp → χχ + j is loop-induced. We show three exemplary Feynman diagrams in Fig. 7. Further diagrams arise by alternatively attaching the final state gluon to another t-, t-or g-line or to the gluon vertex in the upper diagrams.
We calculate the corresponding LHC limits as follows. For the implementation of the model we use FeynRules [93,94] utilizing FeynArts [95] and NloCT [96] to calculate the relevant UV/R 2 counterterms [97]. We generate parton-level events with Mad-Graph5_aMC@NLO [54,98] using the NNPDF 2.1 set [99]. In this context we make use of the loop-induced mode [100] of MG5aMC, which we interface with Ninja [101,102], Golem95 [103] and CutTools [104] for the internal tensor reduction. To gain statistics we apply the parton-level cut p jet T > 200 GeV. We simulate the succeeding parton shower with Pythia 8 [69]. The detector simulation is performed within Check-MATE 2 [105,106] using Delphes [107] where jets are defined via the anti-k T algorithm [108] within Fast-Jet [109,110]. We confront the simulated events with the latest 13 TeV monojet analysis implemented in CheckMATE based on 3.2 fb −1 of data collected by the ATLAS detector [111].
Since the relevant process pp → χχ + j involves at least three heavy particles in the loop, the corresponding cross-section is highly loop-suppressed. More precisely, we find σ(pp → χχ + j) < 10 −5 pb for λ χ = 1, which is seven orders of magnitude smaller than the leading SM background. Therefore, we find that the monojet limits are relevant only for very large values of λ χ , i.e. λ χ > ∼ 7 (for small ∆m) for which the perturbative calculation is already highly questionable. For ∆m = m t the limit is pushed to λ χ > ∼ 9 cf. right panel of Fig. 6 (dark blue shaded region denoted by 'LHC loop-ind.'). This limit can only be improved modestly with new data. For illustration we show the projected sensitivity for 3 ab −1 at 13 TeV where we furthermore optimized the cuts by using TMVA [112] to perform a boosted decision tree analysis [113], providing an estimate for the maximal sensitivity at the LHC (see dark blue, dashed line in the right panel of Fig. 6). Note that the sensitivity, however, does not improve significantly beyond an integrated luminosity of 100 fb −1 due to systematic uncertainties.
For DM masses below m h /2 the invisible Higgs decay h → χχ is open and constitutes another relevant search channel at the LHC. These searches have been performed by the ATLAS [114,115] and CMS [116] collaborations. Here we adopt the 95% C.L. limit BR inv < 0.24 [116] based on an integrated luminosities of 5.1, 19.7, and 2.3 fb −1 at center-of-mass energies of 7, 8, and 13 TeV, respectively. We compute the invisible decay width using the loop-induced hχχ coupling discussed in Sec. III A and use Γ SM = 4.03 MeV [117] to compute BR inv = (1 − Γ SM /Γ inv ) −1 . We do not take into account a possible interference with the dark matter production via direct t/t-loop considered above as we expect the Higgs exchange contribution to dominate for an on-shell Higgs. 5 Furthermore, the selection criteria for Higgs invisible decay searches are expected to further reduce the contribution from the direct t/t-loop. The resulting constraint on the thermal relic scenario is shown in Fig. 8, it excludes a large range of ∆m for DM masses below 53 GeV. The exclusion region is also superimposed in the right panel of Fig. 6 (see blue shaded region labeled by 'Higgs inv.').

E. Searches for long-lived particles
For mediator decay lengths that are comparable to or larger than the size of the LHC detectors the mediator traverses significant parts or all of the detectors. Due to its strong interaction with the detector material the mediator is expected to hadronize and form Rhadrons [118]. At the LHC R-hadron searches are performed exploiting highly ionizing tracks and an anomalous time-of-flight [119][120][121][122]. We use R-hadron searches to constrain both the region of conversion-driven freezeout and the WIMP region. In the entire former region the mediator decay length is large compared to the size of the LHC detector rending the mediator to be detectorstable. The respective limits from the CMS analyses using 18.8 fb −1 of data at 8 TeV [119] and 12.9 fb −1 of data at 13 TeV [120] (preliminary analysis) are shown in Fig. 9 as the dark and light blue shaded area, respectively, and exclude a large fraction of this parameter space. Note that the allowed parameter space in the conversion-driven freeze-out region (after imposing limits from BBN, see Sec. IV F) does not extend above DM masses of about 1.6 TeV. The stop production cross section in the corner of maximal m χ is around 0.11 fb at the 13 TeV LHC [123]. On the basis of the signal efficiencies and background predictions reported in [120] and assuming that the number of observed background events follows its expectation, we conclude that the entire conversion-driven freeze-out region can be probed with an integrated luminosity of approximately 300 fb −1 at 13 TeV.
In the WIMP region the mediator decay-length is typically smaller than the detector size. However, due to the high sensitivity to R-hadron signatures, the respective searches can also impose constraints on intermediate lifetimes for which only a certain fraction of R-hadrons traverse the tracker. Here we use the reinterpretation of the above searches for finite lifetimes presented in [30]. We take the result for the 'generic model' [124,125] and display the 8 TeV limit only (in the relevant region of small masses the limits at 13 TeV are not stronger). Even in the WIMP region R-hadrons probe a small part of the parameter space with small mass splittings close to the boundary of the conversion-driven freeze-out region that is otherwise not robustly constrained by other searches.

F. BBN constraints
The presence of a long-lived, (color-)charged mediator during BBN can affect the predictions for the primordial abundances of light elements in two ways: through energy release from its decay [126][127][128] and through the formation of bound states with the baryonic matter [129,130]. In the present case of a hadronically decaying mediator the former effect provides the stronger constraints due to strong hadro-dissociation processes. In order to es-timate the constraints 6 from BBN we apply the results from [127] for a hadronic branching ratio of 1 presented in terms of the abundance and life-time of the late decaying particle. We use Y t evaluated at x = 0.01 GeV/m χ which is always well after freeze-out (cf. left panel of Fig. 3) and correct for a possible fraction of mediators that have already decayed such that we obtain Y t before its decay. 7 The respective abundances range between Y t = 10 −14 and 10 −13 . We take into account the slight dependence on the mediator mass by linearly interpolating (and extrapolating) the limits provided in [127] for 100 GeV and 1 TeV in log-log space. The resulting constraint on the parameter space is shown in Fig. 9 (red shaded region). It excludes small mass splittings up to around 14 GeV corresponding to mediator lifetimes around 10 s.

V. CONCLUSION
In this work we presented a comprehensive phenomenological study of a simplified DM model where a neutral Majorana fermion is responsible for the observed DM abundance and interacts via a top-philic tchannel scalar mediator. We find that this setup comprises a complex but well-defined phenomenology, giving rise to a large amount and distinct combination of signatures, some of which go beyond those for typical WIMP searches.
The cosmologically viable parameter space encompasses two distinct regions in parameter space in which DM is produced by two different mechanisms: A "WIMP region" where DM freeze-out occurs via DM annihilation or coannihilation, and a region for small mass splitting and DM mass below about 2 TeV where the DM abundance is set by the mechanism of conversion-driven freeze-out. This region in parameter space has not been considered before for this model, and is characterized by a relatively weak coupling of DM to the SM, such that conversions between DM and the mediator are not strong enough to maintain chemical equilibrium. In this case extended Boltzmann equations have to be solved, taking conversion processes into account. For a top-philic mediator we find that the conversion rate is dominated by scatterings while (inverse) decays are suppressed due to the large top mass. For smaller DM masses also 2 → 2 scatterings involving top quarks and W ± bosons become kinematically disfavored, and 2 → 3 scatterings play an important role. Nevertheless, even outside of the "WIMP region" the coupling strength of DM to the SM is still strong enough to thermalize the DM candidate at high temperatures, wiping out any dependence on the initial abundance at the end of the reheating process. Therefore the predictivity of the WIMP paradigm is preserved also in the region governed by conversion-driven freeze-out. In addition we find a peculiar feature of the "relic density surface" in the three-dimensional parameter space allowing for multiple values of the coupling for given DM and mediator masses, due to an interplay of several competing effects.
For the "WIMP region" we have paid special attention to the regime of light DM masses below the top mass, for which the leading annihilation channel to a pair of top quarks is kinematically forbidden. In this case 2 → 3 annihilation processes as well as loop-induced couplings to the Higgs boson and to gluons play an important role. We analyzed direct detection constraints from Xenon1T via loop diagrams, γ-rays observed by Fermi-LAT from dwarf galaxies and the galactic center, antiproton constraints from AMS-02, and collider searches for a colored top-partner as well as invisible Higgs decays, which cover complementary parts of the parameter space and are subject to different types of uncertainties. In addition, we derived limits from loop-induced monojet signatures, that are, however, limited to large couplings. For mass splitting ∆m = m t LHC stop searches exclude the range m χ = 56−400 GeV, and Higgs invisible decay excludes m χ < ∼ 53 GeV. Direct detection is more sensitive for m χ < ∼ 150 GeV and mediator masses m t > ∼ 200 GeV, and in addition closes a small gap between LEP and LHC bounds. For ∆m = m t LZ will probe DM masses up to 2 TeV.
Within the region governed by conversion-driven freeze-out the colored mediator can only decay via 4body processes, and is stable on detector time-scales. This gives rise to R-hadron signatures, which rule out a significant part of the respective parameter space. For a mass splitting between the mediator and DM below ∼ 14 GeV, the lifetime of the mediator exceeds limits from Big Bang nucleosynthesis, leading to complementary constraints and imposing an upper bound on the DM mass of around 1.6 TeV within this region. We estimated that the entire conversion-driven freeze-out region can be probed with R-hadron searches at the 13 TeV LHC with approximately 300 fb −1 of data. Due to kinematic suppression of the mediator decay for light DM masses also a small part of the "WIMP region" can be constrained by R-hadron searches.
While most of the parameter space with light dark matter mass, in particular below the top mass, is excluded, a number of blind spots remain. Apart from the Higgs resonance region, and a small patch for ∆m = m t , the case of almost degenerate masses is notoriously difficult to exclude. To probe this region a dedicated analysis searching for heavy charged and colored particles with non-prompt decay should be considered in the future.