Coloured coannihilations: Dark matter phenomenology meets non-relativistic EFTs

We investigate the phenomenology of a simplified model with a Majorana fermion as dark matter candidate which interacts with Standard Model quarks via a colour-charged coannihilation partner. Recently it has been realized that non-perturbative dynamics, including the Sommerfeld effect, bound state formation/dissociation and thermal corrections, play an important role in coannihilations with coloured mediators. This calls for a careful analysis of thermal freeze-out and a new look at the experimental signatures expected for a thermal relic. We employ a state of the art calculation of the relic density which makes use of a non-relativistic effective theory framework and calculate the effective annihilation rates by solving a plasma-modified Schrodinger equation. We determine the cosmologically preferred parameter space and confront it with current experimental limits and future prospects for dark matter detection.


Introduction
Astrophysical and cosmological observations clearly show that the Standard Model of particle physics (SM) describes only about 20% of the matter content of the Universe. The remaining ≈ 80% can be described to excellent precision by a new cold and collisionless kind of matter. Understanding the nature and the origin of this dark matter (DM) is one of the most pressing issues in high energy physics and cosmology.
Weakly interacting massive particles (WIMPs) produced by thermal freeze-out in the early Universe are arguably one of the most appealing candidates for DM. In this picture the DM has been in thermal equilibrium with the SM plasma at high temperatures. While the Universe cools down, the annihilations rate of DM pairs into SM particles cannot keep up with the expansion rate, and DM finally decouples from the thermal environment. However, the great improvements of direct detection experiments, which exclude cross sections as low as 4.1×10 −47 cm 2 [1], lead to an increasing tension between experimental results and theoretical expectations in models in which the direct detection and annihilation cross sections are related by a crossing symmetry [2][3][4].
There exist a number of well-known scenarios in which this crossing symmetry is not realized. Consequently, the correlation between the annihilation rate and the direct detection rate is more complicated than naive estimates might suggest. One particularly interesting possibility is that the DM is not the only state in the dark sector which is present during freezeout. In this case "coannihilation", i.e. processes which include additional degrees of freedom from the dark sector in the initial or final state, play a role in setting the relic density [5,6]. The thermally averaged (co)annihilation rates which control the evolution of the DM in the early Universe get suppressed very rapidly once the mass difference between the DM and the coannihilation partner is comparable to the temperature T at freeze-out. Therefore, only coannihilation partners with a comparatively large interaction rate or a very strong mass degeneracy with the DM allow for efficient coannihilations. This makes colour charged mediators, which possess a large QCD annihilation cross section, popular coannihilation partners for the DM [7][8][9][10][11][12][13][14].
The phenomenology of coloured mediators is extremely rich. On the one hand, it is well established that non-perturbative effects ranging from Sommerfeld-enhancement to boundstate formation play a non-negligible role in DM production for coloured coannihilations [7,9,[15][16][17][18][19]. As freeze-out occurs in the thermal bath of the hot and early Universe, a careful consideration of the total annihilation rate which takes non-perturbative effects in a thermal background into account is crucial for a reliable determination of the relic density. In particular, the inclusion of bound-state effects for DM annihilation has recently been shown to have a large impact on the relic abundance and increases the largest DM mass compatible with the observed DM density considerably [17,[20][21][22]. On the other hand, coloured mediators can be tested at the LHC, where their QCD interactions allow for copious production, and direct detection experiments. In this paper we try to connect theses different aspects and combine state-of-the-art predictions for thermal freeze-out with a detailed study of the phenomenology at colliders and direct detection experiments. In this way we can confront the cosmological expectations for a thermal relic with other observables and map out the parts of the theoretical favoured parameter space which are in agreement with laboratory searches.
This article is organized as follows. First, we introduce a simplified model for coannihilation which we use for our calculations. In Sec. 3 we discuss the computation of the relic density with a particular focus on the challenges which arise due the presence of a colour charged coannihilation partner in the plasma in the early Universe. Next, we discuss experimental observables which allow to test the model under consideration with present and future detectors in Sec. 4. We map the results of laboratory experiments on the theoretically favoured parameter space of a thermal relic in Sec. 5. Finally, conclusions and outlook are offered in Sec. 6. Some technical details regarding the non-relativistic operators accounting for velocity suppressed annihilation cross section and the thermal potentials used in the analysis are given in the Appendix A and B respectively.

Simplified Model
The simplified model that we consider consists of a gauge singlet Majorana fermion (χ) and a scalar field (η), the latter is a singlet under SU L (2) but carries non-trivial QCD and hypercharge quantum numbers. In the MSSM framework, the Majorana fermion can be identified with a bino-like neutralino and the scalar with a right-handed stop or more generally any right-handed squark. However, we do not fix couplings to their MSSM values and treat them as free parameters.
The Lagrangian for this extension of the Standard Model can be expressed as [23] where H is the Standard Model Higgs doublet, M η the mass of the mediator and M χ the mass of the DM particle. The Yukawa coupling between η and χ is denoted by y while λ 2 and λ 3 are the self-coupling of the coloured scalar and its coupling to the Higgs, respectively. The coupling λ 1 is left for the Standard Model Higgs self interaction and P R (P L ) is the right-handed (left-handed) projector.

Deriving the relic density
With the precise observation of the Cosmic Microwave Background (CMB) by the Planck satellite the cosmological abundance of DM has been measured to be Ω DM h 2 = 0.1200 (12) [24].
Given that this is by far the most precise measurement of DM, an accurate theoretical prediction of the relic density is key ingredient of any realistic study of DM phenomenology. In this context it is crucial to determine the degrees of freedom which are relevant for the production of the DM in the early Universe. In models with coloured mediators the relative abundance of η and χ is typically proportional to the ratio of their equilibrium densities since the conversion rates between them are much larger than the Hubble rate over a large range of temperature and chemical equilibrium holds during the freeze-out process. 3 Consequently, two regimes can be distinguished: (i) for large values of ∆M = M η − M χ the abundance of coloured mediators in the plasma is highly suppressed during the freeze-out of DM and established methods for the computation of the relic density are readily applicable (ii) for ∆M/M 0.2 the abundance of η is non-negligible and can influence the dynamic during freeze-out.
The freeze-out mechanism sets the abundance of the Majorana DM fermions and coloured scalars. The picture is the standard one for heavy particles annihilating in pairs in a thermal ensemble: starting with an equilibrium distribution, the chemical equilibrium is gradually lost when the temperature drops below the heavy particle mass T M χ and recombination processes become Boltzmann suppressed. Eventually the production and annihilation rates cannot keep up with the expansion rate of the Universe and the heavy particles become so rare that they decouple from the thermal medium. The coloured scalars are non-relativistic objects around the freeze-out and during later stages of the annihilations. Indeed kinetic equilibration makes their velocity v ∼ T /M to be smaller than unity. In the case of coannihilating species close in mass with the actual DM particle, the system of Boltzmann equations can be simplified and the evolution of the whole system can be traced with a single effective Boltzmann equation [5,26,27] The total equilibrium number density, which accounts for both particle species of the dark sector (χ and η), is and the effective annihilation cross section reads The mass splitting ∆M T comprise both the in-vacuum and thermal contributions [21], the latter originated from the interactions with particles in the plasma and N c = 3. The sum in eq. (3.3) extends over all coannihilating particles, namely i = χ, η. The key ingredient is the  [29], RT stands for relativistic theory. Right: Landau damping and gluodissociation processes, i.e. 2 → 2 scattering process with light plasma constituents and 2 → 1 process that involves a thermal gluon. In diagram a), dashed lines stand for coloured scalars and the gray blob for a resummed HTL gluon propagator. In diagram b), the solid line stands for a coloured scalar pair in the singlet configuration whereas the double solid line for a colour octet configuration. thermally averaged annihilation cross section in eq. (3.1) which controls the dynamics in the early Universe. Normally, σ eff v is determined by averaging the in-vacuum cross sections σ ij v over the center-of-mass energies in a thermal plasma and reweighting the results according to eq. (3.3). However, the situation can be more complicated and subtle if the coannihilating species can interact with particles from the thermal bath. In our case, the scalar particles feel QCD strong interactions with quark and gluons and the thermally average cross section exhibits a rather strong temperature dependence. Therefore, the in-vacuum result is not a good approximation [17][18][19][20][21][22]28]. In particular bound-state formation is especially efficient at small temperatures and make late-stage annihilations important.
Upon the standard definition of the yield parameter Y ≡ n/s, where s is the entropy density, and the change of variable from time to z ≡ M χ /T , one can rewrite the Boltzmann equation and obtains where m Pl is the Planck mass, e is the energy density, and c is the heat capacity, for which we use values from ref. [30]. We scan over the model parameters y ∈ [0.1, 2] and λ 3 ∈ [0, 1.5] and for every combination (y, λ 3 ) we impose the condition Ω DM h 2 = 0.1200(12) [24]. This way a relation between the DM mass and mass splitting is obtained to constrain the model.

Coloured-scalar dynamics in a thermal bath
The coloured scalars are moving slowly and they can undergo several interactions before annihilating into light Standard Model particles. For example, multiple exchange of gluons leads to the Sommerfeld effect that increase (decrease) the annihilation rate for an attractive (repulsive) potential experienced by the heavy pair. Moreover, the coloured scalar can potentially form bound states in the medium. As the two particles are close to each other in the bound state, their annihilations become more efficient and reduce the relic density. The Majorana fermion does not interact with gluons at tree level, therefore long-range interaction do not appear in χχ annihilations or χη coannihilation processes.
Treating the dynamics of colour charged particles in a medium (including bound-state dynamics formation and dissociation) is not a trivial task. Even though we are interested in coloured scalars in this work, the problem is quite similar to heavy quarkonium in medium which can be described as two quasi-static colour sources sitting in a thermal ensemble. Profiting from the recent progress in quarkonium physics at finite temperature [31][32][33] we adopt an effective field theory (EFT) approach which is inspired by these results. In this framework, tailored quantum field theories are derived to effectively study the physics at a given energy/momentum scale. Two sets of energy scales play a role. On the one hand, the non-relativistic scale M η , the typical momentum M η v and kinetic energy M η v 2 with the hierarchy M η M η v M η v 2 are relevant. On the other hand, there are thermal scales such as the plasma temperature πT and the Debye mass m D which accounts for the chromoelectricscreening inverse length, or equivalently, a thermal mass acquired by the gluons. 4 A schematic arrangement of the energy scales is shown in figure 1. In our case, we are after the dynamics at the binding energy (and kinetic energy) scale, identified with M η v 2 , and its interplay with the thermal scales.
Based on the finding for heavy quarkonium the following in-medium effects are expected to be relevant for coloured scalars: the Coulomb potential is modified by thermal effects and, more importantly, an imaginary part arises in the potential that corresponds to two different processes, namely the Landau damping and thermal break up from bound to unbound pair triggered by a gluon of the plasma [31,32]. The first process involves 2 → 2 scatterings with light particles from the medium whereas the second one describes the absorption of a thermal gluon. In a non-EFT framework, they are known as gluo-dissociation [34] and dissociation by inelastic parton scattering [35]. The equivalence of the thermal breakup process to gluo-dissociation has been analysed in [36], whereas the relation between the Landau-damping mechanism and the dissociation by inelastic parton scattering has been investigated in ref. [37]. We sketch the two processes in figure 1. So far, only one process at 4 πT stands for the temperature scale where the factor π is a remnant of the Matsubara modes of thermal field theory. For coulombic and near-coulombic states, the velocity is v ∼ αs. So the momentum and binding energy/kinetic energy can be written as Mηαs ad Mηα 2 s for the scale estimation that it more appropriate for a bound state, whereas those in terms of velocities better describe scattering states. a time has been included in the bound-state dynamics and relic density determination even though both processes are known to be active in a thermal environment. In the following analysis, we shall consider both processes in the attempt to attain a more comprehensive description of bound-state effects for this model.
In this work, we adopt a formalism suited for addressing the thermal annihilation of nonrelativistic particles [18,19]. The core of the method is to interpret σ eff v as a chemical equilibration rate Γ chem , which can be defined either on the perturbative or non-perturbative level within linear response theory [38]. In this language Γ chem can be traced back to the thermal expectation values of four-particle operators that describe heavy particle annihilations in a non-relativistic EFT (NREFT). In the following, we review the approach with a focus on the model under consideration and the connection with NREFTs and potential NREFTs (pNREFTs).

Generalized Sommerfeld factors from NREFT and spectral functions
The early Universe sets the stage for heavy-particle annihilations. The freeze-out dynamics starts off when the temperature of the thermal bath drops the DM mass. The Majorana fermions and coloured scalars can be considered heavy because their mass is larger than any other energy scale in the problem. In particular, their typical momentum M χ(η) v and kinetic energy M χ(η) v 2 are smaller than the mass. The same stands for the plasma temperature in the regime of interest. In other words, momentum modes as large as the DM mass are highly virtual and irrelevant for the non-relativistic dynamics in a medium with M χ T . Hence, the degrees of freedom are thermalized Standard Model particles with momentum of order T and non-relativistic η and χ particles (they both have similar masses in our setting). In order to obtain the corresponding EFT, one has to integrate out hard energy/momentum modes of order M χ by setting all the other scales to zero, including the thermal scales. 5 As a result, the matching between the fundamental theory (2.1) and the non-relativistic theory can be done as in-vacuum. The so-obtained low-energy theory describes non-relativistic Majorana fermions and coloured scalars interacting with light degrees of freedom. As far as the coloured particles are concerned, the prototype for such EFT is non-relativistic QCD (NRQCD) describing heavy quarks [39]. An EFT describing non-relativistic Majorana fermions has been formulated in refs. [40,41]. As typical of non-relativistic field theories, the number of heavy particles is exactly conserved. Since light degrees of freedom with momenta of the order of the large mass M χ have been integrated out, it seems we cannot describe accurately heavyparticle annihilations into final state particles with hard momenta. However, the inclusive annihilation rate can be recast in terms of an amplitude that conserves the number of the heavy particles thanks to the optical theorem [39,42]. Therefore, particle-antiparticle annihi-lations are accounted for by introducing four-particle operators. The operators are organized in a 1/M 2 χ expansion, where the small parameter is the velocity of the heavy particles. Such an expansion resembles the velocity expansion of the annihilating states (s-wave, p-wave and so on). The matching coefficients of these operators contain the effects of the relativistic degrees of freedom and are obtained by matching four-point Green's functions between the fundamental and the low-energy theory, here the NREFT.
Starting from the relativistic theory written in (2.1) and expanding in the non-relativistic limit the scalar and fermion fields, the Lagrangian that comprise the four-particle dimension-6 operators read [21] L d=6 where φ, ϕ and ψ annihilate a scalar particle, a scalar antiparticle and a DM fermion respectively; α, β are colour indices in the fundamental representation of SU(3) and p, q = 1, 2 spinorial indices. The optical theorem relates the matching coefficients, which fix the parameters of this first low-energy theory, to heavy scalar annihilations into light degrees of freedom, i.e. the in-vacuum cross sections. Neglecting the masses of SM particles they read [21] where C F = (N 2 c − 1)/(2N c ) and h stands for the Yukawa coupling between the SM Higgs and the quarks. In this model, the Majorana DM pairs annihilate in two quarks and it features a helicity-suppressed s-wave rate [23,43,44]. This is reflected by c 1 = 0 in eq. 3.6 because the masses of the Standard model particles are set to zero. However, in order to treat the boundary between the coannihilation regime and the standard freeze-out correctly an accurate description of χχ annihilations is highly desirable. Therefore, we will include the leading contribution to χχ annihilations even though it is of higher order in 1/M χ . The situation is different for heavy and light quarks respectively. In the former case, we compute the matching coefficient for the first four-particle operator in (3.5) keeping a nonvanishing top-quark mass, and we find where m t is the topquark mass. In the latter situation, the ratio m u(d) /M χ is tiny and we take the p-wave contribution as the leading term for χχ annihilations. In the NREFT language, we have to include 1/M 4 χ suppressed operators in (3.5) which give velocity suppressed contributions to the in-vacuum cross section of order v 2 /M 2 χ . When thermally averaging these terms, a T /M suppression appears [27]. The dimension-8 operators and their matching coefficients are given in Appendix A. We stress that we do not include the effects of either the top-quark mass or p-wave suppressed contributions for χ-η coannihilations and coloured scalar pair annihilation as these already possess a non-vanishing contribution at lowest order in 1/M χ .
As shown in refs. [18,19] the thermally averaged annihilation cross section can be then written as where · · · T stands for the thermal average and O i are the four-particle operators listed in eq. (3.5). This expression makes the factorization of the scales M χ M χ v, M χ v 2 , πT, m D manifest. On a physical ground, hard annihilations happen at a distance scale much smaller than any other introduced by the soft energy scales. These processes are local and unaffected by the medium. Soft momenta are taken into account in the thermal average of the very same operators and characterize the dynamics of the coloured scalar before their annihilations: multiple gluon exchange, mass corrections, colour-phase decoherence. A linear response analysis shows that [18,19] where E is the energy of relative motion of the pair after factoring out the center of mass dynamics, and Λ is a cutoff restricting the average to the non-relativistic regime, in particular . 6 The main advantage of this approach is that the spectral functions ρ i contain all the dynamical information of the annihilating heavy pair [18,19]. As far as the coloured scalar are concerned, a repeated gluon exchange leads to the Sommerfeld effect for the scattering states and also to the formation of bound states for the colour-singlet channel. At finite temperature, both possibilities need to be considered simultaneously, and thermally averaged over the whole phase space. After reshuffling some terms we can define thermally averaged Sommerfeld factors [18,21,45] which read as follows: where N i are the number of contractions for each operator. Note that theS i factors account for all non-perturbative effects in-medium and encode Sommerfeld enhancement/suppression, bound state formation and thermal effects. In terms of the latter quantities, the annihilation cross sections can be written as 7 where the thermally averaged Sommerfeld factors for processes with χ in the initial state are unity since the DM is a singlet under QCD. A possible way to determine the spectral functions for the non-relativistic pairs (and then theS i through eq. (3.8)) is to solve a plasma modified Schrödinger equation. It implements a resummation of repeated interactions between the coloured scalars and the particles from the medium. The spectral function has been shown to correspond to the imaginary part of the Green's function that solves the very same Schrödinger equation [46,47] where the Hamiltonian has the standard form is an in-medium potential whereas Γ T (r) represents the interaction rates, i.e. real scatterings with the plasma constituents (Landau damping and gluon emission/absorption) that has been interpreted as the imaginary part of the heavy-pair potential [31,32]. We remark that the potentials are functions of the temperature and their explicit form depends on where the thermal scales sit with respect to the non-relativistic scales (see figure 1). The potential and the width in eq. (3.11) comprise an r-independent and r-dependent part. The former comes from self-energy diagrams of the each particle of the pair, whereas the latter comes from a gluon attached to the two different heavy particles. The heavy scalar potentials in the different annihilation channels can be interpreted in a rigorous and modern language as matching coefficients of a pNREFT that we discuss in the next section.

QCD-potentials for coloured scalars in pNREFT
In this section we deal with the non-relativistic heavy-scalar pairs and the QCD potentials affecting them. The quantum field theory that describes the energy modes of order M η v 2 is a pNREFT. We are interested in such energy modes since we want to study the boundstate formation and dynamics. Well-established examples of effective theory of this kind are pNRQCD for a heavy quark-antiquark pair at zero temperature [48,49] and its generalization to finite temperature [32] (for the Abelian version see refs. [50,51]). In practice one has to integrate out the soft momenta of order M η v characterizing the gluon-momentum transfer between the heavy scalar-antiscalar pair. These modes are indeed still dynamical in the NREFT. 8 Then, depending on whether the temperature scale is larger or smaller than M η v, different pNREFTs are obtained (see refs. [32,33,52] for heavy quarkonium). In our case, the degrees of freedom are the singlet, octet and sextet wave function fields and the matching coefficients correspond to the singlet, octet and sextet potentials respectively. The presence of the Majorana DM fermion allow for the annihilation channel ηη → qq, so that a particleparticle potential is also relevant.
As mentioned in the previous section, the details of the low-energy theories depend on the assumed hierarchy of scales. In an expanding Universe, the temperature of the plasma is decreasing with time and different arrangements of the energy scales are expected to occur. For example, the condition πT M η v, which is valid around the early stages of the scalar annihilations is not going to hold at later stages and one should plug different potentials in the Schrödinger equation (3.10). As far as the potential arising from inelastic parton scattering, results are available for quarkonium in the case M η v πT [32,37], however one would also need the interpolating regime M η v ≈ πT that is realized during the early Universe cooling. The real part of the potential in the latter situation has not been calculated yet and is subject of a future work [53]. Here, we do not pursue a sophisticated analysis by inserting many different thermal potentials in eq. (3.10) because the focus of this work is on the model phenomenology and the ties with the relic density constraints that capture bound-state formation. In particular, we use the high temperature potentials derived in ref. [21] that rely on the Hard thermal loop (HTL) approximation and account for a Debyescreened Yukawa potential and soft 2 → 2 scatterings with particles of the plasma (i.e. the Landau damping as shown in figure 1, diagram a). We exploit these potentials for the whole temperature range. Even though this is not strictly rigorous, the corresponding effects are decreasing with the temperature and capture the physics of the thermal environment in the QCD sector. As a refinement with respect to the previous study undertaken in ref. [21], we add the gluodissociation contribution for the singlet channel in the thermal width Γ T (r) (process b in figure 1) and the accompanying thermal modification in the real part of the potential V T (r). We use the expressions derived in ref. [32] in the static limit. In the pNREFT language, this process is induced by chromolectric transitions between singlet and octet fields. It is expected to be subleading at temperatures πT, m D ∆V , where ∆V is the difference between the singlet and octet energy. Additional details on the pNREFT and the thermal potentials are offered in Appendix B. In this way, we include both the processes relevant for bound state dynamics in a unified framework.
Let us now write down the effective Lagrangian. In the high-temperature assumption, namely πT M v, the temperature scale has to be integrated out first. The main outcome is that the gluon propagator is modified by a thermal mass m D , and also an imaginary part. The latter accounts for real scatterings with the parton constituents of the plasma (quarks and gluons). The gluon propagator takes the so-called Hard Thermal Loop (HTL) form [54,55]. Then the typical heavy-scalar momentum transfer can be integrated out and we obtain a pNREFT HTL (we use the same language as given in refs. [56,57]). The degrees of freedom of such EFT are colour singlet, octet and sextet field. The leading contribution to the Lagrangian in the static limit reads where the trace refers to a colour trace, then S, O and Σ are the singlet, octet and sextet fields, are the r-dependent potentials whereas δM s(o,Σ) are the self-energy corrections for the heavy pair and L HTL is the HTL Lagrangian for the gluons and light quarks. These quantities contain both a real and an imaginary part and read [21,31,32] where the auxiliary function Φ is given by [31] Φ (3.17) The real and imaginary part of eqs. (3.13-3.16) stand for the r-dependent and r-independent terms of V T (r) and Γ T (r) and have to be inserted in the plasma-modified Schrödinger equation (3.10) for each annihilation channel. In contrast to the case of heavy quark-quark pairs studied in [58], the antitriplet field does not contribute in the case of scalar pairs due to Bose-Einstein statistics. We would like to stress that even though the potentials have previously been derived in [21] their interpretation in terms of the pNREFT (3.12) for heavy coloured scalar comes as a novel aspect in this work. By working with a pNREFT we can make contact with the formalism developed for heavy quarkonium at finite temperature which (i) facilitates the identification of the relevant degrees of freedom at the energy scale M η v 2 (ii) allows a rigorous derivation of the thermal potentials according to the assumed hierarchy of scales (iii) provides a unified framework in which all the processes relevant for bound-state dynamics are included.

Qualitative discussion of relic density
Let us briefly discuss our results for the generalized Sommerfeld factors and the corresponding effect on heavy scalar annihilations. In the left panel of figure 2, we show the generalized Sommerfeld factorS 3 with and without the effect of gluo-dissociation/bound-state formation by radiative emission and thermal corrections to the real part of the potential for M η = 3 TeV (see eq.(B.1) for details). At high temperatures we find a very small effect on the Sommerfeld factorS 3 that is slightly reduced by the additional dissociation process. On the contrary, at smaller temperatures that corresponds to z > ∼ 300, the contribution from O → S + g (octet into gluon and singlet) is sizeable and it helps in forming bound states via the process of the radiative emission as described in other works [17,20,22]. 9 In this work we focus onS 3 that enters the thermally averaged cross section in the singlet channel. Here, the bound-state formation boosts the annihilation rate at small temperatures and gives a large effect. For the repulsive channels, we use the values ofS 4 andS 5 as derived in ref. [21]. Their effect is to reduce the corresponding octet and sextet annihilation channels, however the effect is less pronounced than in the singlet channel. In our work we do not include the chromoelectric transitions that affect the octet self-energy (there are no transitions of the sextet field because the antitriplet configuration is absent for scalars). Nevertheless this should be pursued for a more rigorous theoretical treatment and it poses an interesting topic for future research [53]. 10 In the right panel of figure 2, the relic density is displayed for y = 0.3 and λ 3 = 0.1 and for selected values of ∆M/M . The dashed lines correspond to a computation where only the Debye screened Yukawa potential and the thermal width by the inelastic parton scattering is 9 In these works the spectral function method is not exploited. Instead, the bound-state formation cross sections calculated by the authors of refs. [17,20,22] are related via the Milne relation to the dissociation (or ionization) cross sections which are obtained from the dissociation rates ΓT . 10 The octet self-energy in the hierarchy of scales MQv T mQv 2 mD, ΛQCD has been computed in ref. [59] for heavy quarkonium, where MQ is the heavy quark mass.
included [21], whereas the results denoted by dotted-dashed lines also account for the 2 → 1 process for bound-state dissociation/formation. We notice that the effect is more prominent as the mass splitting decreases. Indeed, the final relic abundance is more sensitive to the coloured scalar annihilations when the scalar particles are very close in mass to the DM fermion. For comparison we also show the in-vacuum results for vanishing mass splitting. As can be seen the DM abundance is severely overestimated even for ∆M/M χ = 0 if the corrections are not taken into account.
Besides the gluon exchange, the Standard Model Higgs boson can induce an attractive potential between heavy coloured scalars in this model. The Higgs exchange has been recently considered in refs. [60,61] and the corresponding effect depends on how the interaction between the coloured scalar and Higgs boson is realized. On the one hand, for a trilinear vertex hη † η induced after electroweak symmetry breaking by the Higgs vev (as in the Lagrangian given in 2.1) the Higgs exchange plays little role [61]. The main reason is the heavy mass suppression that, together with the scalar coupling values considered here λ 3 ≤ 1.5, generates an effective coupling α eff ≡ λ 3 v h /(8πM ) ≈ 0.01 that is typical of weak interactions strengths for the masses under consideration. Therefore, we neglect Sommerfeld enhancement and bound-state formation induced by the Standard Model Higgs boson in the present work. On the other hand, a trilinear interaction which is proportional to a new mass scale can have a sizeable impact if the new scale/interaction strength is large enough [60].

Experimental constraints
Naturally, one wonders whether the cosmologically preferred regions of the parameter space can be probed experimentally. In particular direct detection experiments and collider searches are known to have a great sensitivity to models with colour-charged mediators [62][63][64][65][66] 11 In contrast to freeze-out, which is largely insensitive to the flavour of the quark with which the DM interacts, the rates and signatures at colliders and direct detection experiments depend quite sensitively on the specific quark. In order to bracket the range of possibilities we consider two representative choices in the following: DM interactions with (i) the up quark as a representative of the valence quarks and (ii) the top quark.

Collider searches
In collider searches the impact of the quark flavour on the phenomenology is twofold. First, new production modes which are sensitive to the couplings and masses of the DM can con- 11 Current indirect searches using Fermi-LAT data [67] or antiprotons [68,69] are not sensitive to the canonical annihilation cross section expected for a thermal relic unless astrophysical enhancement factors boost their sensitivity. Since the annihilation rate of coannihilating DM in the Universe today is generically suppressed these searches do not constrain the theoretically favoured parameter space and we will not discuss them in detail.
tribute significantly to the total production rate at the LHC if the DM couples to quarks with a large parton luminosity. Second, the flavour of the quark which is produced in the decay of the mediator η has a non-negligible impact on the sensitivity of experimental searches.

Valence quarks
The production of new physics at colliders receives contributions from different processes. First, gluon-gluon and qq initial states allow for the production of ηη † pairs by gluon exchange due to the colour charge of the mediator. The production amplitude from quarks-antiquark pairs also includes a term with t-channel χ exchange which can dominate over the QCD production for large values of y. In addition, ηη (η † η † ) and ηχ (η † χ) final states are produced exclusively by new physics. Same sign scalar pair production, which is possible because of the Majorana nature of χ, is mediated by t-channel χ exchange. In the regime where new physics contributions to scalar pair production are relevant the ηη rate is typically enhanced compared to ηη † and η † η † due to the large u-quark parton distribution function (PDF) and dominates over the other channels for large y [66]. In addition, (anti-)quark gluon initial states allow for the production of η ( †) χ final states. The amplitude of this production mode is linear in the new physics coupling y and therefore most important at intermediate values of y since the t-channel processes scale with a higher power of y.
At colliders the final states of models with colour charged mediators consist of jets from the decay of the mediator and pairs of DM particles. Since the DM can not be observed with the LHC detectors χ production leads to a momentum imbalance of the observed particles and the signal consists of missing transverse energy (M ET ) in association with jets. In the coannihilation regime the mass different ∆M is small and the jets from η decay are typically soft which makes them hard to distinguish from QCD backgrounds. Consequently, hard partons in the final state are necessary to provide additional handles on the event and suppress the backgrounds. In this case the DM will recoil against the visible jets (or other SM particles) and thus increase the missing transverse energy. Numerous searches for these types of signatures exist, [70][71][72][73][74][75] with monojet searches typically being the most sensitive.
We implement the simplified model in FeynRules [76] and simulate the collider signal with MadGraph5 aMC@NLO [77]. The results are then passed to Pythia8 [78] for showering and hadronization. Detector effects are taken into account with Delphes 3 [79] and we use the Checkmate 2 [80] implementation of the ATLAS monojet search [70] 12 to derive limits on the cosmologically favoured parameter space.
Our analysis shows that the pure QCD production cross section of η mass-degenerate with χ is excluded for M η < 320 GeV. In the mass degenerate regime the new physics production modes increase the reach of the LHC search to 540 GeV for y = 1 and 1100 GeV for y = 2.
However, the relic density constraint imposes a non-negligible mass splitting between the DM and η for these couplings. As a consequence, the region of parameter space which is tested by the ATLAS search does not overlap with the region of parameter space considered in our analysis of the relic density.

Top Quarks
For η interacting with top quarks the story is rather different. First, the parton luminosity of top quarks is completely negligible at the relevant center of mass energies. Consequently, only the pure QCD production of the mediator is relevant in this case. The decay of the mediator is more complicated than in the case of light quarks. In the region of parameter space for which coannihilations are relevant m t + M χ > M η such that the two-body decay η → χt is kinematically not accessible. Therefore, three-body and four-body decays into a b-quark, an on-or off-shell W boson and χ dominate the width and the finals states become rather complex compared to η with couplings to light quarks. The production and decay modes of η are similar to the simplified topologies used in searches for supersymmetry [70,[82][83][84] and the results for stops apply to the model considered here. Currently, the LHC is only sensitive to M χ 500 GeV and the regions of parameter space of interest remains unconstrained.

Direct detection
Direct detection experiments probe the scattering rate of DM particle off nuclei in lowbackground detectors. In the last years the sensitivity of direct searches has seen a huge improvement and the best current limits probe scattering cross sections as low as 4.1 × 10 −47 cm 2 [1]. As a consequence these experiments are now able to test scenarios which feature a suppressed scattering rate such as coannihilations.
In the model under consideration here, the direct detection cross section is rather sensitive to the fundamental interactions between the DM and the quarks. On the one hand, DM interacting with light quarks couples directly to the valence quarks in the nucleus and a sizeable interaction is already generated at tree-level. On the other hand, the tree-level interactions between DM coupling to heavy quarks and the nucleons are highly suppressed since the expectation value of heavy quark pairs in protons and neutrons is negligible. Consequently loop-induced couplings between the DM and light quarks and gluons become relevant and dominate the DM nucleon cross section in this case. In general both spin-dependent (SD) and spin-independent (SI) interactions are generated. However, due to the higher experimental sensitivity SI DM-nucleon scattering dominates the constraints and we have checked that SD scattering is not competitive. Therefore, we will focus on SI interactions in the following.
In general, the SI direct detection cross section on nucleons is given by [85] σ To make the discussion more transparent we will therefore discuss the two cases separately. Before discussing direct detection in detail it should be noted, however, that the limits and prospects of direct searches depend on details of the astrophysics of DM which are not fully understood. In our work we assume the standard halo model, i.e. a truncated Maxwellian velocity distribution, with v 0 = 220km/s, v esc = 544km/s and v e = 232 km/s and taken the local DM density to be ρ 0 = 0.3 GeV/cm 3 in agreement with the standard choice made by the experimental collaborations for the presentation of direct detection limits [1].

Valence quarks
If the DM interacts directly with light quarks the dominant contribution to the effective DM nucleon coupling arises from the tree-level exchange of η, see Fig. 3 for a representative diagram. In addition, the loop-induced coupling between the DM and the Higgs can become relevant for large λ 3 and this lead to an additional contribution to the DM nucleon coupling. Taking all contributions into account one finds [23,86,87] where m h is the SM Higgs mass. The f N T i are related to nucleonic matrix elements of the quark and gluon operators by f N T q = N |m qq q|N /m N and f T G = 1 − q=u,d,s f N T q and parametrize the contributions of quarks and gluons to the nucleon mass. The quantities q(2) andq(2) denote the second moments of the quark and anti-quark parton distribution functions, respectively. In our numerical analysis we use the default values for f T q from micrOMEGAs 5.0 The first term in Eq. 4.2 arises from an effective four fermion interaction χχqq between the DM and the quarks. For DM coupling to light quarks this interaction is generated by tree-level η exchange with a strength given by [86] g q = − y 2 8 In addition, triangle diagram with η and the quarks running in the loop induced an effective DM Higgs coupling g hχχ . Higgs exchange contributes to the χχqq interaction and induces a coupling to the gluons once heavy quarks are integrated out. Due to the small Higgs coupling to valence quarks only the diagrams with the Higgs coupling to η contribute appreciable to g hχχ and the interaction with the quarks can be neglected. In the zero-momentum transfer limit relevant for direct detection the effective coupling is given by where x η = M 2 η /M 2 χ and the quark masses have been neglected. This contribution is always subdominant for DM but can amount to a O(10%) correction in certain regions of parameter space if λ 3 is large.
The last term in eq. 4.2 is also generated by tree-level η exchange and the strength of the interactions is controlled by the same parameter g q which enters in the first term.

Top quark
For DM interacting with heavy quarks, the direct tree-level coupling to the nucleus is absent since the abundance of heavy quarks in the nucleons is highly suppressed. Consequently, higher order processes become important. Two terms contribute substantially to the effective DM nucleon interaction: (i) the effective DM Higgs coupling leads to an interaction via Higgs exchange; (ii) box loops with quarks and η in the loop generate an effective coupling between the DM and the gluon content of the nucleons, see Fig. 4 for representative diagrams contributing to the loop induced coupling.
In the low energy limit, the ensuing DM nucleons coupling reads [23,86,87] where G(2) is the second moment of the gluon parton distribution while b and g G parametrize the strength of addition contributions to the DM nucleon interaction which were negligible for valence quarks. The first two terms arise from the effective DM Higgs coupling and are similar to the ones introduced in the previous section. Due to the large Yukawa coupling of heavy quarks, additional diagrams with the Higgs coupling to the internal quarks line contribute to g hχχ . Taking the full quark mass dependence into account the effective coupling in the zero-momentum transfer limit reads [9] g hχχ = 3y 2 8π 2 where The remaining terms arise from box-diagrams which generate an effective coupling between the DM and the gluon content of the nucleus. The third term is induced by the χχG a µν G aµν operator where G a µν denotes the field strength tensor of QCD. The strength of this interaction is given by [86] (4.9) The I i (M η , m q , M χ ) are comparatively lengthy loop functions which are given in the Appendix of [88]. The box diagram also contributes to an operator which connects the DM to the gluon twist-2 operator, i.e. the traceless part of the energy momentum tensor ΔM/M χ Figure 5: The left (right) panel summarized constraints on a thermally produced DM coupling to u-quarks via a coloured mediator η for λ 3 = 0.0 (λ 3 = 1.5). The red region is excluded by results of the Xenon1T experiment. The orange region shows the part of parameter space which can be tested by a future Darwin-like direct detection experiment. In the upper gray area y ≥ 2 is required to reproduce the correct relic abundance. In the lower gray region the thermal abundance is controlled exclusively the gauge interactions and the observed relic abundance can not be achieved for DM from freeze-out. In the brown region, the binding energy of the lightest ηη † bound state is smaller than 2∆M . The dashed blue lines indicate the required value of y while the dotted black line indicates results obained with the in-vacuum cross section and y = 0.3.
The coefficient of the effective Lagrangian is given by [86] (4.11)

Results
In this section we combine our results for the abundance of thermal relics with the limits from direct detection experiments and LHC searches. Since the relic abundance is known very precisely we can fix one of the model parameters (y) in terms of the others (M χ , M η , λ 3 ) by imposing Ω DM h 2 = 0.1200(12) [24]. In order to make the discussion more transparent we consider different slices of the parameter space and discuss their phenomenology. To be more specific, we choose representative values for λ 3 , i.e. λ 3 = 0 and λ 3 = 1.5, which bracket the impact of λ 3 and allows for an intuitive visualization of the results. In our analysis we focus on ∆M/M χ ≤ 0.2 since coannihilations are subdominant at larger mass splitting and we restrict ourself to M χ ≥ 500 GeV. For lower M χ the late stage annihilations due to boundsstate formation make the freeze-out process sensitive to temperatures at which the QCD potentials become non-perturbative. In this regime our results are no longer reliable and one would need to determine the thermal expectation values on a lattice. Such an analysis goes beyond the scope of this work. As in the previous section, we first discuss the phenomenology of light quarks before turning to top quarks. In order to avoid repeating similar results we limit ourself to DM coupling to up quarks in the following and do not discuss the other light quarks in detail.
As can be seen in Fig. 5 the parameter space for successful freeze-out is limited from above and below. The blue dashed curves reproduce the observed relic density in the (M χ , ∆M/M χ ) plane for representative values of y ranging from y = 0.1 to y = 2.0 in steps of 0.1. They are obtained via the numerical solution of the Boltzmann equation for the DM yield given in eq. (3.4). The effective annihilation cross section σ eff v is defined in eq. (3.9) and it includes the dynamics of coloured mediators in medium, i.e. Sommerfeld enhancement, bound-state formation/dissociation and thermal effects as explained in Sec. 3.2 and 3.3. On the one side, as ∆M/M χ increases coannihilations become less efficient and the value of y which reproduced Ω DM h 2 = 0.1200(12) increases until the limit is reached. On the other side, for low ∆M/M χ coannihilations become very efficient and since some of the coannihilation rates are exclusively set by the Standard Model charges of the mediator we find a region (lower gray area) in which the coannihilation rate is too large too accommodate Ω DM h 2 = 0.1200(12) irrespectively of y 13 . As can be seen the slope of the relic density curves changes noticeably at high and at low ∆M/M χ . The first change at ∆M/M χ in the range 0.1 − 0.2 corresponds to the onset of non-negligible coannihilation and marks to upper edge of the parameter space in which the in-medium effects considered here are relevant. The change of slope at ∆M/M ≈ 5 × 10 −3 in Fig. 5 and Fig. 6 is due to very efficient bound-state formation that corresponds to large coloured scalar annihilation rate at small temperatures. This was anticipated in Fig. 2 where we show the relic abundance for different values of ∆M/M . The smaller the mass splitting the longer the coloured scalars are as abundant as the DM and, therefore, the abundance of the overall dark sector is driven by large annihilation rates of the coloured particles. For comparison, the back dotted line indicates the exemplary results obtained with a pure invacuum computation that neglects Sommerfeld and bound state effects for y = 0.3. By increasing λ 3 to a value close to the perturbative limit (we show results for λ 3 = 1.5) the coannihilation processes get markedly more efficient and the region of freeze-out DM shifts upwards in the M χ , ∆M/M χ plane. This is due to the fact that the Higgs contributes to the singlet channel annihilation rate which is enhanced by the especially largeS 3 Sommerfeld factor in the thermally averaged cross section (3.9). In the brown region, the binding energy of the lowest-lying bound states is smaller than two times the in-vacuum mass splitting, i.e. 2∆M < |E 1 |, such that the lightest two-particle states in the dark sector are the bound states formed by the coloured scalars. Due to efficient chemical equilibration rates that convert η particles into DM particles, it seems plausible that almost all DM fermions convert into the scalars and subsequently annihilate away, so that the model is most likely not able to explain the viable observed DM abundance in this parameter region [21].
Collider searches are not able to probe a part of the cosmologically favoured part of the parameter space for coannihilating DM. The ATLAS monojet search is not sensitive enough the exclude the pure QCD-monojet cross section, i.e. η pair production in association with a jet, for m η ≥ 320 GeV. The additional production modes due to new physics which depend on y enhance the monojet cross section substantially and dominate the monojet rate for large y but even in this case the cosmologically preferred η masses for a given y are outside of the reach of the LHC search by several hundred GeV.
Direct detection experiments can test the coannihilation scenario and, for λ 3 = 0, a big part the parameter space is already excluded by the null-results reported by XENON1T [1]. In addition, the prospects for testing this scenarios with future experiments are excellent and a detector with an exposure similar to the proposed DARWIN experiment [89] could probe the entire region. However, the larger mass splitting associated with a larger value of λ 3 decrease the DM nucleon cross section and makes these scenarios harder to access experimentally. Nevertheless, even for the case of λ 3 = 1.5 current experiments are already starting to exclude parts of the high and low mass region and not even the tip of the coannihilation strip, i.e. M χ ≥ 10 TeV, remains beyond the reach of future experiments.
In the case of DM interacting with top quarks the picture is somewhat different. On the one hand, the cosmological preferred parameter space is essentially insensitive to the quark flavour and remains unchanged. On the other hand, the collider and the direct detection phenomenology is quite different. As detailed in Sec. 4.1 collider searches are currently not sensitive to M χ ≥ 500 GeV irrespectively of the other model parameters and, therefore, the LHC limits do not appears in Fig. 6. Also direct detection experiments currently lack the sensitivity necessary to probe the parameter space of thermal DM coupling to tops. However, this situation will improve substantially once future experiments begin to collect data. A DARWIN-like detector can test a significant part of the parameter space. In contrast to the case DM coupling to up quarks the direct detection prospects for DM coupling to tops become better as λ 3 increases and for λ 3 = 1.5 it is almost possible to reach M χ = 10 TeV for large values of y. This is due to the fact that the loop mediated DM Higgs coupling, which is sensitive to λ 3 , contributes significantly to the direct detection rate while the suppression due to the higher mass splitting associated with large λ 3 is less pronounced than for valence quarks. All considered the prospects for thermal DM produced by coannihilations with a colour charged mediator in the future are encouraging and strengthen the case for an ultimate direct detection experiment such as DARWIN.

Conclusions and Outlook
In this paper we connect the recent improvements for the relic density determination of WIMPs with current experimental limits and future prospects. To be specific, we consider a simplified model with a Majorana fermion as DM candidate and a strongly interacting scalar. The latter acts as a mediator between the Standard Model and the dark sector and is assumed to be close in mass with the DM in order to allow for significant coannihilations. In this case, the DM relic density is influenced by the annihilations of the coloured scalars that feel strong interactions. Due to repeated gluon exchanges and interactions with light plasma constituents in the thermal medium, the determination of the coloured-scalar annihilation rates is nontrivial and requires care. In the last years a great effort has been made to include the effects of Sommerfeld enhancement and bound-state formation in the relic density computations. These new results indicate that previous analyses underestimated the impact of coannihilation on the relic density and call for a reassessment of the experimental capabilities and a detailed study of the parameter space preferred by thermal freeze-out. It is manifestly subtle and complicated to include bound-state dynamics in a thermal medium due to an intricate interplay between non-relativistic and thermal energy scales. In this work, we adopt a formalism that connects the thermally averaged cross section with the determination of the spectral function for non-relativistic particle pairs. The latter can be extracted by solving a plasma-modified Schrödinger equation. In doing so, one can describe the two-particle state entering the hard annihilation and account for the dynamical formation of bound states in a thermal bath, without any assumption on the nature of the annihilating pair. As far as the coloured scalar pairs are concerned, different processes for bound-state dissociation/formation are active: 2 → 2 soft scattering with plasma constitutes and gluodissociation/absorption. These are captured by the imaginary part of the thermal potentials experienced by the coloured scalars. In addition, the real part of the potential is modified by the thermal plasma. Profiting from recent progress in the context of heavy quarkonia in a hot QCD medium, we use a pNREFT for heavy coloured scalars in the high-temperature limit. Here, the thermal potentials for the singlet, octet and sextet fields are understood as matching coefficients, and account for Debye screening and inelastic parton scattering in the real and imaginary part respectively.
In order to attain a more complete description of the coloured scalar dynamics and make contact with other results in the literature, we have also taken into account the process of bound-state dynamics triggered by a thermal gluon. In the pNREFT language, it is described by singlet to octet transitions. When extracting the generalized Sommerfeld factors we find results which agree with the general arguments provided by the EFTs for heavy quarks in a medium and the corresponding power counting. Singlet-to-octet transitions give a substantial contribution for small temperatures (here identified with z > ∼ 300) whereas they plays little role at higher temperatures. To the best of our knowledge, both the processes responsible for bound-state dynamics have not been considered simultaneously in this model before.
The thermal potentials among non-relativistic particles are a key ingredient to a better understanding of DM dynamics in a thermal bath. However, it is also important to focus on the rate equation governing the time evolution of the DM abundance in order to connect it with phenomenology. Here, we used a Boltzmann equation and plug in it a thermally average cross section that tracks the nature of the two particles annihilating: scattering or bound states depending on the temperature and the other model parameters.
With these ingredients a precise determination of the DM relic abundance is possible and the cosmologically favoured regions of the parameter space can be mapped out. We focus on regions with strong coannihilations, i.e. ∆M/M ≤ 0.2, and restrict ourself to M χ ≥ 500 GeV. An analysis of M χ < 500 GeV would be of considerable interest, however, for such low DM masses QCD can enter the non-perturbative regime before the freeze-out process is complete. The formalism that we use can handle the non-perturbative regime of QCD, however one would have to evaluate the thermal expectation values of the heavy particles on a lattice. The preferred regions of parameter space can then be confronted with direct detection and collider experiments. Unfortunately, we find that current LHC searches are not sensitive to the region of parameter space of interest here due to the small y/high M η required to reproduce the observed relic density. In contrast to collider searches, direct detection experiments have an excellent sensitivity to coannihilating DM. In models of DM interacting with light quarks current bounds on the DM nucleon scattering cross section already exclude significant parts of the parameter space and future detectors such as the proposed DARWIN experiments will be able to probe the complete parameter space conclusively. The situation is a bit less optimistic if DM interacts preferably with heavy quarks and current experiments are not yet able to exclude DM interactions with top quarks. Nevertheless, even in this case future experiments will be able to probe a significant fraction of the parameter space.
Finally, we would like to stress that there remains room for improvements in the determination of the thermal relic abundance. First of all, a more detailed classification of the thermal potentials according to different assumption on the scale hierarchies would be highly desirable and the impact of changing scale arrangements during the evolution of the system ought to be explored on a quantitative level. Moreover, the dynamics induced by chromoelectric transitions for the colour octet configuration have not been included so far and an assessment of the corresponding effect on the repulsive thermal potential could improve the precision of the relic density calculation further. Recently, an alternative derivation of the rate equations for the DM number density has been carried out [90]. In this case, an ab initio out-of-equilibrium treatment is pursued and a useful systematization of the different rates governing the DM dynamics can be handled. The general equations can be simplified according to different assumptions and, for example, the result obtained from a linear response theory is recovered as a limiting regime close to chemical equilibrium. Perhaps, another interesting and promising approach to a detailed dynamics of color-singlet and -octet densities is the one discussed in refs. [59,91,92]. Here the heavy quark-antiquark pair is interpreted as an open quantum system in a thermal environment and an out-of-equilibrium dynamics is also taken as a starting point. It would be worth exploring this direction for applications to relic density determination for coloured coannihilators and compare quantitatively the outcomes of the different theoretical approaches.
terms in the Lagrangian are where the included operators read The arguments of the 2S+1 L J indicate the angular momentum state of the χχ pair which is created or annihilated by the operator, ↔ ∂ is the difference between the derivative acting on the spinor to the right and on the spinor to the left, and T (ij) stands for the symmetric traceless component of a tensor, i.e. T (ij) = (T ij +T ji )/2−T kk δ ij /3, that appears in eq. (A.5) through the derivative and Pauli matrices vector components. At variance with NRQCD there is no covariant derivative because the DM Majorana fermion is a singlet under QCD. Moreover, the number of independent operators is smaller than NRQCD because of the relation σ k pq σ k rs = 2δ ps δ qr −δ pq δ rs as noticed in ref. [21]. Another difference is that a Majorana fermion pair is even under charge conjugation and following the discussion in ref. [93] one can already select which operators are allowed or not by looking at the corresponding J P C .
The matching is done by computing the one-loop diagrams in the fundamental theory given in figure 7. One has to expand and keep terms up to order (v/M χ ) 2 and match onto the corresponding amplitude in the NREFT. As mentioned in the body of the paper, we only need the imaginary part of the one-loop amplitude. We perform the calculation in the center of mass frame for the incoming χχ pair. The incoming Majorana fermions have momenta p and −p, whereas the outgoing states have p and −p . By energy conservation we have |p| = |p | and the velocities are defined as v = p/E and v = p /E, with E = p 2 + M 2 χ . The result reads (we drop the spinors and we use the notation as in ref. [39]) This quantity has to be matched with therefore we obtain In this calculation the quark mass is set to zero. The number of contractions of the dimension-8 operators contain two powers in the velocity. When computing the corresponding thermally averaged cross section, one obtains a thermal integral that appears in the velocity expansion (i.e. non-relativistic expansion) of the cross section [27]. Finally it reads [21,23] This is the thermally averaged cross section that we use to obtain the relic density constraints in the case of DM interacting with light quarks.

B EFTs and thermal potentials
Here we give more details about the numerical implementation of the potentials to extract the spectral functions and the corresponding generalized Sommerfeld factors. As mentioned in the body of the paper, we do not pursue a detailed analysis of the various arrangement for the non-relativistic and thermal scales. We restrict to a simplified scheme with the aim to include the two different processes that describe bound-state formation/dissociation and the corresponding effects on DM annihilations. We highlight the similarities and differences with respect to the numerical derivation and setting adopted in ref. [21]. First, for the numerical evaluation of the Schrödinger equation one needs α s for a wide range of energies and the Debye mass for a large range of temperatures. We remind the reader to the appendix of ref. [21] for details on how these quantities are fixed since we stick to the same implementation. We just notice that even going down to temperatures smaller than 1 GeV, we can rely on a perturbative determination of the potentials. Our observation is based on a recent lattice study [94] that suggests the absence of string-tension contributions for temperatures T > 2 T c , where T c ≈ 150 MeV is the QCDcrossover temperature. Since we consider the smallest DM matter mass to be M χ = 0.5 TeV and we integrate the Boltzmann equation down to z = 10 3 , we are within such a case.
Second, as mentioned in the body of the paper, we add the contribution of the gluodissociation process in V (r) T and Γ(r) T for a better estimation of the thermally averaged Sommerfeld factors. In the pNREFT language, this amounts at considering singlet to octet transitions mediated by a chromoelectric gluon. We take the following expression as calculated in ref. [32] that are valid for the scalar field η in the fundamental and antifundamental representation of SU(3) c Here, P stands for the principal value, n B is the Bose-Einstein distribution and ∆V is the difference between the singlet and octet static potentials ∆V ≈ N c α s /(2r). For a numerical estimate of ∆V , we used the inverse Bohr radius for the inverse distance, 1/r = 1/a 0 = M α s C F /2 and then ∆V ≈ M α 2 s . In eq. (B.1) the typical expansion of pNRQCD appears, namely the small distance r [32,48,49]. When solving the Schroedinger-like equation (3.11) with the addition of the terms in eq. (B.1), we impose them to contribute when 1/r ≥ πT is satisfied. This is needed to justify the usage of that potential in first place and also to avoid spurious effects in the spectral function determination at large distances. Finally, we notice that the Coulomb potential at small distances (and temperatures) is already recovered in the HTL expressions by summing the r-dependent with the r-independent part in eqs. (3.13) and (3.16). Despite it is not rigorous to use the HTL potentials at small temperatures, namely πT < 1/r, they still provide a well-defined behaviour of the solution at both large and small distances.