Impact of Sommerfeld Effect and Bound State Formation in Simplified $t$-Channel Dark Matter Models

The existence of a dark matter model with a rich dark sector could be the reason why WIMP dark matter has evaded its detection so far. For instance, colored co-annihilation naturally leads to the prediction of heavier dark matter masses. Importantly, in such a scenario the Sommerfeld effect and bound state formation must be considered in order to accurately predict the relic abundance. Based on the example of the currently widely studied $t$-channel simplified model with a colored mediator, we demonstrate the importance of considering these non-perturbative effects for correctly inferring the viable model parameters. We emphasize that a flat correction factor on the relic abundance is not sufficient in this context. Moreover, we find that parameter space thought to be excluded by direct detection experiments and LHC searches remains still viable. Additionally, we illustrate that long-lived particle searches and bound-state searches at the LHC can play a crucial role in probing such a model. We demonstrate how future direct detection experiments will be able to close almost all of the remaining windows for freeze-out production, making it a highly testable scenario.


Introduction
One of the prominent outstanding puzzles of particle physics and cosmology is the origin and nature of the most abundant matter component of the universe, namely Dark Matter (DM). Apart from its manifestation through gravitational effects affecting several astrophysical observations of "ordinary" matter (assuming the DM has a particle physics origin), we do not have a clear picture about its nature. For decades, one of the most promising candidates for DM have been 'weakly interacting massive particles' (WIMPs), cold non-relativistic particles that are able to interact with the primordial thermal bath of Standard Model (SM) particles (and potentially with other species beyond the SM) in the early universe. In this period of time, they can annihilate into SM particles and, as the universe expands and cools, they freeze-out leaving behind a relic that we can observe today as a clear imprint on, for example, the spectrum of the cosmic microwave background [1]. A number of beyond the standard model (BSM) physics scenarios predict particles that can serve the purpose of WIMP dark matter, allowing us to formulate strategies to try to detect them [2].
Over the last few decades, a large variety of experiments have been searching for Dark Matter (DM), ranging from the production of dark sector particles at colliders [3] to direct detection of DM-nucleon scattering on earth [4,5] as well as the observation of astrophysical signals of annihilation and/or decay products of dark sector particles [6]. A conclusive positive (direct) signal for DM is still missing, while the data from these experiments can be used to constrain models of DM.
One possibility for why WIMP dark matter has evaded detection so far, could be the existence of a rich dark sector. Such scenarios open up for the possibility of colored coannihilation [7][8][9][10][11][12][13][14][15], which naturally lead to a heavier expected dark matter mass, evading current experimental searches. These types of models are currently intensively studied by the theoretical and experimental community in form of so-called t-channel simplified models [16][17][18]. In particular, the interplay of the aforementioned experimental searches and the requirement not to overproduce dark matter allows for complementary constraints on the model parameters.
If the DM-SM coupling is sizable, the evolution of DM proceeds via the thermal freezeout of DM, which typically takes place at temperatures of T ∼ m DM /30 when DM itself is non-relativistic. In scenarios where dark sector particles interact via a light mediator, long-range effects such as the Sommerfeld Effect (SE) [19,20] will enhance (diminish) the annihilation cross-section [21] due to the resulting attractive (repulsive) potential, altering the theoretically predicted dark matter abundance [9,[22][23][24][25].
Additionally, dark sector bound states can radiatively form via the emission of the light mediator particle. Their decays effectively deplete the dark-sector particle densities, acting as an additional effective annihilation channel for DM and therefore affecting the theoretical prediction of the dark matter relic density [26][27][28][29][30][31][32]. Based on the methods developed in [33], it has been shown that Sommerfeld enhancement and bound state formation (BSF) arising from non-Abelian interactions can lead to O(50% − 300%) corrections to the theoretically predicted dark matter abundance [34][35][36]. Given this sizeable correction, we demonstrate that taking into account the aforementioned non-perturbative effects is an imperative task in order to derive correct exclusion limits in colored simplified t-channel models. We also stress that the same argument could be extended to any model featuring colored co-annihilation or similar conditions.
In particular, we show, compared to previous works on this subject, that several important effects which had not been previously taken into account significantly expands the scope of both the model phenomenology and the complementary ways of experimental detection. Specifically, we point out the following salient features of this work: • So far, previous work considering non-perturbative effects in co-annihilation scenarios assumed that the DM annihilation cross-section is solely determined by the annihilations of the mediators charged under a non-Abelian gauge group [34,35]. We also take into account the effects of a non-zero coupling of the mediator to the DM candidate, which is mandatory for annihilations of the mediator to impact the DM relic abundance. Furthermore, the presence of this coupling is crucial to test the model at various experiments. We provide a detailed analysis of the interplay of current and future experimental searches and point out the areas of parameter space where bound-state effects are large.
• When computing the BSF rate, we include the effects of the three-gauge-boson vertex, as described in [34], which was neglected previously [27,35]. Additionally, we estimate the effects of a non-negligible mediator decay width on the efficiency of DM depletion induced by bound-state effects.
• In [37], two of the authors of this work analyzed the simplified t-channel models described in this paper by calculating direct detection and collider limits at next-to-leading order (NLO) and showed that spin-independent direct detection limits at NLO can be significantly more important than the corresponding tree-level spin-dependent part. However, in the aforementioned work, the co-annihilation scenario and the impact of the Sommerfeld enhancement as well as of bound state formation were not taken into account in the calculation of the DM relic density. The present paper aims at showing the full impact of both these improvements.
• While recently the effects of BSF on the relic density have been considered in the context of the conversion-driven freeze-out [38] and the superWIMP mechanism [39], we investigate in this article, the scenario where DM is in thermal equilibrium with the dark sector and the SM.
• In order to derive exclusion limits, we fix the DM-SM coupling by the requirement not to exceed the observed dark matter abundance. This differs from previous works on the subject [16,17], where the DM-SM coupling was fixed requiring the mediator decay width to equal 5% of its mass. The latter approach, however, would overestimate the experimental conclusion limits in the mass-compressed region, which is the main focus of this work.
For this purpose, we focus on a class of simplified t-channel DM models called u R , d R and q L [37] (or S3M uR, S3M dR and S3M QL, in the spirit of [16,17]), where DM is a Majorana fermion χ and has a Yukawa interaction via a color triplet scalar X with a SM quark. Interestingly, the model parameters, the coupling of DM to Standard Model (SM) particles g DM , the DM mass m DM and the mass splitting in the dark sector ∆m, are constrained by a combination of cosmological, collider and direct detection limits. In this work, we demonstrate how the corresponding experimental exclusion limits and expected experimental signatures are altered when considering the effect of Sommerfeld enhancement and bound state formation on the DM relic abundance.
We find that: • When including SE, we find that the maximally allowed DM mass and mass splitting in the co-annihilating region of the parameter space are shifted from (m DM , ∆m) (1 TeV, 30 GeV) to (m DM , ∆m) (1.4 TeV, 40 GeV), leading to a correction around (40%, 33%).
• BSF provides a significant correction to the annihilation cross-section when g DM < g s .
In this regime, it adds up to the SE, leading to a stronger impact and must be included in order to obtain a legitimate result. Since BSF effectively acts as an additional annihilation channel, it always increases the effective annihilation cross-section. We find that including the SE and BSF leads to a shift of the maximally allowed DM mass and mass splitting from (m DM , ∆m) (1 TeV, 30 GeV) to (m DM , ∆m) (2.4 TeV, 50 GeV), hence a correction around (140%, 66%).
• Depending on the parameters considered, corrections from the SE can be positive or negative, while BSF always increases the annihilation cross-section. As a result, a simple flat correction factor ("K-factor") is not sufficient for an approximate consideration of these effects.
• In case of an observation, considering SE and BSF is important in order to infer the correct model parameters (m DM , ∆m, g DM ).
• SE and BSF extend the parameter region that would lead to underabundant thermal dark matter, where the observed DM relic density can only be produced by non-thermal DM production. Consequently, this region might not be fully probed by long-lived particle searches in contrast to the expectation without the inclusion of non-perturbative effects.
• We find that bound state searches at the LHC offer a unique opportunity to constrain couplings in the region of 10 −6 g DM 10 −2 , which otherwise typically evades the constraints from prompt collider searches, direct detection or long-lived-particle searches.
The analysis leading to these conclusions is organized as follows: In section 2, we give a detailed description of the simplified t-channel models analyzed and the relevant processes for the relic abundance calculation. In section 3, we illustrate some theoretical aspects of SE and BSF in the simplified t-channel model and we discuss their impact on the relic density and the model parameters in 4. In section 5, we summarize the constraints utilized from spin-independent and spin-dependent searches, while in section 6, we explain how we exploit prompt collider searches, including the search for BSF at the LHC, and long-lived particle signatures. Finally, in section 7, we present our combined results and we elaborate on the interplay of the various constraints and their potential to exclude parts of the parameter space. Most importantly, we discuss the impact of SE and BSF on the estimation of the correct exclusion limits. Moreover, we show the corresponding projected exclusion constraints from future experiments and highlight the potential reach of long-lived particle searches and of searches for dark sector bound states at the colliders. We conclude in section 8.

Simplified t-channel models and Dark Matter Cosmology
In this section, we briefly describe the t-channel simplified model [16-18, 37, 40] as well as the various thermally averaged cross-sections that are relevant for evaluating the DM relic abundance. The t-channel model we consider consists, in addition to the SM, a SM-singlet Majorana fermion χ which is the lightest dark sector particle, and three color-triplet complex scalar fields X i (i indicates the generation) which interact with χ and the SM quarks via a Yukawa coupling g DM . The scalars are charged under the SM gauge group (SU (3) × SU (2)) Y and its simplest form there are three possible quantum number assignments possible: The three possible choices of the mediator's quantum numbers correspond to three different models, which we label as the u R , d R and q L models, respectively. The dark sector features a Z 2 symmetry such that χ is the lightest stable particle and our DM candidate. The interaction Lagrangian of the dark sector particles is thus given by: where D µ is the covariant derivative and the index i runs over the quark and mediator flavours of the model considered (up-type right-handed quarks, down-type right-handed quarks and left-handed quarks). P L and P R are the left and right handed projectors respectively. The Yukawa couplings,g DM , are chosen to be real valued, flavour-diagonal and flavour-universal for simplicity, implying that g DM,ij = g DM δ ij 1 . processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. processes contributing to (co-)annihilations described in Tab. ??. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model.

1
(i) Co-annihilation Figure 1: Representative Feynman diagrams corresponding to the main subset of processes contributing to (co-)annihilations described in Tab. 1. For simplicity, we don't show here possible interfering diagrams from crossing symmetries (for example, we would have an u-channel for χχ → qq, XX † → qq and XX → qq). We also only illustrated the gluon gauge vertices, since the strong coupling dominates. In diagram 1h, we could also have the interaction of any other SM gauge boson with the quarks, while in diagrams 1b, 1c, 1d, 1e and 1i, the gluon can be replaced by a photon or a Z boson in the u R and d R models and additionally by W ± bosons in the q L model. As a result of its Yukawa interaction 2 , the DM number density depends both on direct pair annihilation of DM, particles, χχ → SM SM, as well as co-annihilation and colored annihilations processes into SM particles involving χ − X, χ − X † , X − X † and X − X as initial scattering states. The latter determining the density of the scalar mediators. A representative class of Feynman diagrams are shown in Fig. 1.
Under the assumption that all Z 2 -odd particles will finally decay into dark matter and will be in equilibrium with each other until freeze-out [46], we track the evolution of the total dark-sector comoving number density (or yield )Ỹ =ñ/s, whereñ is the total dark-sector number density and s is the entropy density of the universe, resulting in the sum of the comoving yields of the co-annihilating Z 2 -odd species: We can then write an effective Boltzmann equation forỸ as a function of the variable x = m DM /T in the following form: where c = π/45 m Pl m χ , (2.5) with g χ = 2 and g X = 3 being the internal degrees of freedom of the Majorana particle χ and the colored scalars X, m Pl the Planck mass, g * (g * S ) the number of effective relativistic degrees of freedom for the energy (entropy) density of the Universe and Hereby, the effective annihilation cross-section in Eq. (2.4) is given by where σ ij v ij comprises of all the annihilation cross-sections of two co-annihilating species i and j. If the scalars are much heavier than DM, their abundance gets quickly Boltzmannsuppressed and the only relevant process for determining the DM density is the direct χ − χ annihilation. On the other hand, when δ 1, the scalar mediators continue interacting for a longer time and their thermally-averaged (co-)annihilation cross-section contributions to Eq. (2.10) are significant even around the freeze-out of the DM particle candidates. We stress that in order to derive Eq. (2.10), two crucial assumptions are made. First, the rate of elastic scatterings of dark sector particles with the SM bath is much larger than their annihilation rates, due to the fact that the number density of the SM particles are not Boltzmann-suppressed. This allows us to approximate f i ∼ f eq i , where f i is the phase space density of the particle species i. In addition to this, DM and its co-annihilating partners are considered to be in chemical equilibrium, meaning that the processes that convert a dark sector particle into another one are faster than the expansion of the universe. In other words, the inter-conversion rate χ ↔ X is larger than the Hubble rate, which allows us to approximate Y i /Ỹ Y eq i /Ỹ eq 3 . For concreteness, Eq. (2.10) can be explicitly written by means of the allowed co-annihilation processes, enlisted in Tab. 1 and depicted in Fig. 1, as where the multiplicities in front of the various terms come from the sum in Eq. (2.10) and by considering that scalars and antiscalars give the same contribution (we will omit the † in processes other than for XX † pair annihilations). From this expression, we can appreciate how, for δ 0, colored annihilation and co-annihilation processes can provide a sizable contribution to the total annihilation cross-section. On the contrary, for large δ, these processes are negligible so that the t-channel pair annihilations of DM particles into SM quarks is the dominant channel. In the degenerate limit, δ = 0, the relative weights of Eq. (2.11) respectively amount to (2.12) In the u R and d R models, we consider three copies of SU(2)-singlet scalars, each of which carrying only colored degrees of freedom; this implies that g X = 3 for each scalar, while for the Majorana singlet we only have two spin configurations and g χ = 2. The q L model, instead, features SU(2)-doublet mediators; this implies the existence of six copies of colortriplet scalars which make the denominator in Eq. (2.12) even larger, resulting in a smaller effective annihilation cross-section. In the co-annihilating region, we will see that this will be the main difference between the q L model and the u R and d R models in computing the relic abundance (cf. Sec. 4). While the relative contributions to the total effective crosssection depend on the internal degrees of freedom of the particles involved, the individual Process Contribution to σv v rel Color Structure BSF Table 1: Summary of the most relevant processes contributing to the (co-)annihilation of DM (the corresponding Feynman diagrams are depicted in Fig. 1). In the second column, we present the coupling structure combined with the suppression factor caused by the mass splitting δ in the dark sector. Here, g gauge generically represent the gauge coupling in the vertex qqA, with A being the involved gauge boson (these are g, γ, Z, and W ± , although the latter applies only to the q L model). In the XX † → qq process, α and β generically indicate the contribution from the Yukawa-and the gluon-mediated process. The third column shows the dominant velocity dependence of σv rel in the limit of low relative velocities. A distinction between massive and massless quark is also highlighted, with the latter being the most relevant case for our analysis. The fourth column gives the decomposition of the squared matrix element into the individual contributions of different color configurations [R], with f 1 and f 8 being some functions of the coupling constants. The fifth column indicates the ability of the initial-state particles to form () or not to form () a bound state and a corresponds to the suppressed BSF when these are not in conjugate representations.
cross-sections can have very different magnitudes.
In Tab. 1, we highlight the contributions to σ eff v rel as a function of the couplings involved and of the relative mass-splitting δ (second column). Moreover, we show the dominant velocity-dependent contribution in the limit of small relative velocities (third column) 4 . In the co-annihilating parameter region, we see that the t-channel direct annihilation of DM into quarks as well as the pair annihilation of the scalar-antiscalar triplets into quark-antiquark pairs are velocity-suppressed in the limit of massless quarks and therefore subdominant for almost all the parameter space we are interested in. Clearly, for lighter DM and small masssplittings, the top quark mass would not be negligible. Nevertheless, as we will see (cf. Section 4), typical Yukawa couplings in this region of the parameter space are g DM 10 −1 and the process is therefore suppressed with powers of g 4 DM . In fact, as a general feature, in those regions of the parameter space where the Yukawa coupling g DM is smaller than the strong coupling g s , processes that do not involve vertices with gluons are negligible even if they are not velocity-suppressed (e.g., XX → qq annihilations).
Finally, as we will thoroughly describe in the next section, the individual thermallyaveraged cross-sections can receive large corrections from non-perturbative effects generated by long-range, colored interactions between color-charged, initial states, so that the actual contribution in percentage of one channel can be quite sizable. For this reason, in the fourth and fifth columns we also show the color structure of the amplitude and the possible contributions from BSF to the effective annihilation process, as we will explain in more detail in the following subsection.
3 Long-range effects in the simplified t-channel model Processes that involve two massive, colored particles in the initial states are subject to the long-range behaviour of the gluonic Coulomb potential generated by the strong force exerted between them. This influence distorts their scattering wavefunctions in such a way that they can experience a further attraction or repulsion, depending on the color charge they carry, that is to say on their color representation, which directly influences the gauge coupling α g related to the gluon potential (see the definition later in Eq.3.3). The impact of this effect becomes important when the incoming particles travel at non-relativistic speeds. The typical interaction distance is given by the inverse of the average relative momentum k = µv rel exchanged between two scattering particles (µ being their reduced mass). When the interaction distance approaches the Bohr radius (µα g ) −1 of the interacting system, implying 4 In this work we avoid referring to the term in the velocity expansion of σv rel that has a velocity dependence of v 0 , v 2 , v 4 , ... as (s, p, d, ...)-wave and so on. This terminology is borrowed from the expansion of the probability amplitude into partial-waves with different values of the orbital angular momentum , M = M , where the terms with = (0, 1, 2, ...) are called (s, p, d, ...)-wave and so on. However, it is inconsistent to identify a term in the velocity expansion of σv rel with a certain partial-wave, as it has been discussed in detail in [35]. In fact, the momentum expansion of the matrix element for a given partial-wave results in M ∼ n p +2n .
Consequently, one finds σv rel ∼ |M| 2 = |M | 2 ∼ ,n,n p 2( +n+n ) [35,50]. This implies that terms proportional to p 2 are not purely "p-wave" and actually receive contributions from states with = 0 (e.g., when n = 1 and n = 0). For the successive terms in the momentum expansion scaling with p 2m , the contribution emerges from m + 1 different partial-waves so that utilizing the partial-wave terminology fails for m = 0. This clarification is of particular importance when dealing with Sommerfeld corrections to processes that have no velocity-independent terms [35]. α g ∼ v rel , non-perturbative effects from the continuous exchange of gluons between the two interacting states become relevant. This is the so-called Sommerfeld effect [19]. Moreover, the presence of an attractive potential can naturally lead to the formation of bound states of the color-charged particles that also affect interaction rates. In the following sections we summarize the pivotal results related to these effects in a non-Abelian gauge theory and apply them to our simplified model.

Color decomposition and Sommerfeld effect
Let us indicate with R 1 and R 2 the color representations under which two incoming particles transform. Their interaction can then be decomposed into a direct sum of irreducible representations, The gluonic interaction in the non-relativistic regime is described at leading order by a static Coulomb-like potential that depends on the irreducible representation as follows Here, Q represents the averaged momentum transfer in the interaction. The effective gluon coupling constant in a given irreducible representation of the unbound scattering state α

3)
C 2 (R) being the quadratic Casimir invariant of the given representation R. Notice that the k [R] factor can be either positive or negative, meaning that the potential can be either attractive or repulsive. We evaluate the coupling constant α at the scale of the average transferred momentum of the scattering states, Q = µv rel , by considering corrections to the RGE of the β-function of α s (Q) up to NNLO, as given in [51] and as employed by numerical codes used for our scans (e.g., in micrOMEGAs [52]). For simplicity, we suppress in the following, the dependence on the average transferred momentum Q, while consistently taking it into account in our numerical calculations. Since the scalar fields of our models belong to the fundamental representation 3 of SU(3) (the antiscalar fields to the conjugate one,3), we can easily list the different possible potentials we obtain from Eq. 3.2. Since 3 ⊗3 = 1 ⊕ 8 and 3 ⊗ 3 =3 ⊕ 6 (the case with3 ⊗3 gives identical results to the last one since C 2 (R) = C 2 (R)), we have the following potentials: .

(3.4)
We can see how the singlet state accounts for the most attractive potential. In this work, we neglect the effects of the Yukawa potentials generated by electroweak bosons. For light DM, the Yukawa exponential suppression makes the influence at large distances negligible. In fact, in order for long-range effects to have a sizable impact, the Bohr radius needs to be smaller than the inverse mass of the force mediator, (α m X /2) −1 m −1 A . For typical electroweak couplings of order α ∼ 0.02, effects of the electroweak potential start to be important for m X 8 TeV. However, thermal freeze-out typically takes place at T ∼ m DM /30, which results in T ∼ m X /30 in the coannihilating regime. Since electroweak symmetry breaking (EWSB) occurs at T EWSB ∼ 150 GeV, electroweak gauge bosons are massive during freeze-out only if m X 4.5 TeV. Thus, below this threshold, effects of the electroweak potential can be safely omitted 5 .
If the dominating term in the bare annihilation cross-section σ 0 is velocity-independent, the Sommerfeld effect results in a simple multiplicative factor ( multiplying the "s-wave" term, which is the only one in the partial-wave expansion) [35,55]: is the s-wave Sommerfeld factor and where c [R] is a factor coming from the color decomposition of the amplitude and marks the relative contribution of theR initial state to the total crosssection 6 . For the potentials illustrated in Eq. (3.4), one finds [35] [6] . (3.10) The case of annihilation into a quark-antiquark pair is more involved, as it is not straightforward to decompose the color structure of the amplitude, due to the fact that both QCD-and Yukawa-dominated diagrams contribute and interfere. In general, one needs to compute the single color contribution of each term in the squared matrix element and then determine the color weight factors, that we indicate here as general functions f [1] and f [8] of the couplings involved (e.g., g s and g DM ). Since in our model we encounter this situation only for processes that are velocity suppressed across the vast majority of the parameter space, we choose to neglect the Sommerfeld effects for them. In the last process in Eq. (3.10) we consider the case of distinguishable particles both in the initial and in the final states. This is the situation, for example, for a X i + X j → q i + q j , where i and j refer to two distinct flavours. By employing Eq. (3.6) in our model, the four types of Sommerfeld factors needed are: For the Coulomb potential, the function S 0 (ζ s ) is given by rel and, depending on its sign, results in either an enhancement (S 0, [1] and S 0, [3] ) or a decrease (S 0, [8] and S 0, [6] ) of the perturbative, tree-level cross-section. Therefore, as soon as v rel α g, [R] , the Sommerfeld factor can have a sizable impact on the cross-sections involved, so that the annihilation processes we discussed in Sec. 2 are affected significantly.
For annihilation processes that are dominated by velocity-dependent terms (see Table 1) we should also take into account the Sommerfeld corrections to the different partial-wave contributions. In fact, one can show (see, e.g., [35,55,56]) that each term of the squared amplitude at a given must be multiplied by Therefore, as described in footnote 4, at a given order in the momentum expansion, one should in principle consider for each partial-wave term the corresponding S . For instance, in the v 2 rel term of the expansion of σv rel one needs to consider both an s-wave and a p-wave term, which are corrected by S 0 and by S 1 = S 0 (1 + ζ 2 s ), respectively. In what follows, we will only take into account the dominant s-wave Sommerfeld corrections and will omit these higher partial-wave contributions, since the latter come with factors that are equal to 1 plus a correction ζ 2 s ∝ α 2 s , which is of order O(10 −2 ).

Bound-state formation, ionization and decay
In a similar regime where the Sommerfeld effect is relevant, color charged particles can also form unstable bound states (BS) via the emission of a gluon (cf. Fig. 2): Due to the non-Abelian nature of the strong force, gluon radiation can be emitted both from the external legs (i.e., the wavefunctions of the colored particles forming the bound state), or from a gluon propagator exchanged between the two incoming states, what we call a radiative vertex, see Fig. 3. The wavefunction of the bound state is determined by solving its corresponding Schrödinger equation. This results in a discrete energy spectrum, with energy eigenvalues given by the binding energies E n m = −κ 2 /(2µ n 2 ) with n being the principal quantum number and µ the reduced mass of the two-particle system 7 .
The strong coupling is evaluated at the scale of the average momentum exchange of the particles in the process. In particular, in the ladder diagram (cf. Fig. 2), the average momentum transfer between the bound-state wavefunctions is given by the Bohr momentum . For the strong coupling in the radiative vertex in the BSF diagrams in Fig. 3, α BSF s, [R] , the expectation value of the momentum exchanged is the same as the one of the emitted gluon |P g | = ω, with ω being the gluon energy. From energy-momentum conservation, in the non-relativistic regime, ω must be approximately equal to the difference between the relative kinetic energy of the initial scattering states and the binding energy of the final bound state, neglecting their total kinetic energies 8 .
Finally, for the emission of radiation directly from a mediator exchanged by the scattering states, as shown in Fig. 3, the momentum transfer is Q NA |p − q| 2 . Since k << κ when BSF is relevant, one can approximate this last expression as Q NA κ [R] and α NA s, [1] α B s, [1] . For more details on the computation of the transition amplitudes, we refer the reader to [33,34], on which the formalism employed in this work is based on. Note that the effects of BSF on the relic density have also been considered in the context of non-relativistic effective field theories including non-zero temperature corrections [32,44,[57][58][59][60][61][62].
As described earlier in Sec. 2, the dark-sector (anti)scalars are color (anti)triplets. From group algebra, the relevant decompositions for a scalar anti-scalar interaction into irreducible representations of SU(3) for the radiative transitions are 3 × 3 = 1 + 8, 3 × 3 = 3 + 6 and 8 × 8 = 1 S + 8 S + 8 A + 10 A + 10 S + 27 S , such that the only allowed capture processes for a X ( †) scalar-antiscalar initial state are (X + X † ) [8] → B(XX † ) [1] + g [8] , Here assume that the two-particle initial state is either in the singlet or the octet representation. The final-state particles are the bound state B which belong to an irreducible color representation and the emitted gluon g, is in the adjoint representation 8, by definition. The combination of color algebras in the final states must match the one in the initial state. For example, we see that when the initial state is a color singlet, we can only form an octet bound state. The additional gluon that is emitted, combines with the bound state to form a symmetric singlet representation 1 S . However, as noted earlier, only the singlet potential is attractive, and hence only the first of the processes above in Eq. 3.16 is relevant 9 . For particle-particle interactions the following bound state formation processes are possible (analogous antiparticle-antiparticle interactions exist that we do not reproduce here): (X + X) [3] → B(XX) [6] + g [3] , (3.17a) (X + X) [3] → B(XX) [3] + g [3] , (3.17b) (X + X) [6] → B(XX) [3] + g [6] , (3.17c) (X + X) [6] → B(XX) [6] + g [6] .
In this work, we only take into account the effects of singlet bound states, since they constitue the dominant effects (cf. Eq. 3.4). In fact, the BSF cross-section of a ( [3]) bound state is subject to cancellations arising on the level of the squared matrix element of the formation process caused by the color structure of the process 10 . Hence, in the following, we will focus on the dominant effect of the capture into the color-singlet ground state (n m) = (100) (see Eq. (3.16a)). Generally, excited states, can also open additional annihilation channels through their direct decays into radiation or via bound-to-bound transitions and subsequent decays. We neglect these effects in the present work, and leave them for future improvement 11 . For a scalar-antiscalar pair transforming in the fundamental representation and with degenerate masses m X , the bound-state formation cross-section reads as [34]: Here, α BSF s, [1] represents the strong coupling constant from the gluon emission vertex, while α B s, [1] arises from the ladder diagrams involving the bound state wavefunction. We will omit the subscript [1] in the following. The S BSF function, arises from the overlap integrals involving the scattering-state and the bound-state wavefunctions and reads as Here, ζ S ≡ α S g /v rel and ζ B ≡ α B g /v rel parameterize the ratios between the strong coupling 9 The octet bound state could still have a significant role when a further attractive long-range interaction is present. One example are Higgs-mediated processes in the limit of heavy initial state particles. In fact, being the Higgs a real scalar, it always leads to an attractive force, which counteracts (enhances) the repulsiveness (attractiveness) of the gluonic potential [42,43]. 10 We illustrate the effect of these cancellations in Appendix A. 11 During the completion of this manuscript, two preprints including and discussing excited-state contributions on DM abundance appeared [38,63]. As shown there, the qualitative prediction is not altered significantly, apart for the zero mass-splitting case. For example, in the single-generation dR model analyzed in Ref. [38], the shift in the largest possible DM mass is found to be at the percent level for mass splittings ∆m few GeV, while only for ∆m = 0 the same amounts to O(25%). and the relative velocity for the scattering (S) and bound state (B), respectively 12 . The first parenthesis corresponds to the s-wave Coulomb Sommerfeld factor (cf. Eq. (3.12)), coming from the normalization of the scattering wavefunction, while the second term accounts for its p-wave correction (cf. Eq. (3.13)). The presence of this second factor is due to the fact that the radiative capture is, at leading order, a dipole transition between the initial scattering state and the final bound state. As it is well-known, these transitions impose a selection rule ∆ = ±1 on the orbital angular momentum, such that in order for the final bound state to be in the ground state (100), the initial scattering state must be in a = 1 state, hence the p-wave correction to the Sommerfeld effect.
Bound-state effects give their peculiar contribution in the last factor, which accounts for the convolution of the bound-state wavefunction with the radiative vertex. Similarly to the simple Sommerfeld factor, at large velocities we have S BSF ∼ ζ 4 B ∼ (α B g /v rel ) 4 1 and the cross-section in Eq. (3.18) gets suppressed; at low velocities we obtain again the typical S BSF ∼ α S g /v rel scaling when the interaction between the scattering state is attractive, such that ζ S 1 and ζ B 1. The BSF cross-section in this case is enhanced and can compete with the (co-)annihilation processes. However, when the interaction is repulsive and ζ S 1, the BSF probability becomes exponentially suppressed, S BSF ∝ exp(−2π|ζ S |), as it becomes harder to form bound states when two incoming particles experience a strong repulsive force. Overall, the BSF cross-section is maximized when α S g /α B g = 0.5 and when α ∼ v rel . We may define a thermally-averaged BSF cross-section as where ω = µ/2 (α B g ) 2 + v 2 rel is the energy emitted by the radiated gluon, which equals its kinetic energy minus the (negative) binding energy E 100 = −µ(α B g ) 2 /2 of the newly formed bound state, and f g (ω) = (exp(ω/T ) − 1) −1 is the gluon distribution function, accounting for the Bose-enhancement factor 1 + f g (ω) from the final state gluon. This term is necessary to ensure the detailed balance between bound-state formation and ionization reactions, i.e. ionization equilibrium.
In fact, if bound-states are successfully formed, they can still be ionized by energetic gluons in the thermal plasma and dissociate into their constituents, or can directly decay into radiation. The former process usually dominates at temperatures larger than the binding energy, when particles in the thermal plasma are sufficiently energetic to disrupt the bound state. The ionization rate of a bound state can be written as 12 See also Tab. 1 in Ref. [34] for the various definitions of the strong gauge coupling constants.
where the relation between the ionization cross-section σ ion and the BSF cross-section follows from the Milne relation (see for example Appendix D of [34]). Here, g X = 3, g g = 8 and g B are the internal degrees of freedom of the scalar triplets, the gluons and the bound-states, respectively, and where µ = m X /2 in our model. For the capture into the singlet-state, we simply have g B,[1] = 1. Importantly, we see that the ionization rate becomes exponentially suppressed for T ω, as already anticipated. The decay rate Γ dec,[R] of a = 0 bound state in a given representation R into gauge bosons is computed by taking the s-wave perturbative annihilation cross-section times relative velocity of the corresponding scattering states and multiplying it with the bound state wavefunction evaluated at the origin: Considering the color-singlet ground state decaying into a pair of gluons, one obtains from Eq. (3.18) and from 13 |ψ The formation and subsequent decay of bound states involving dark sector particles can therefore affect the DM relic abundance by effectively opening up a new annihilation channel. In fact, these effects must be incorporated as additional terms into the system of coupled Boltzmann equations that controls the evolution of the number densities of bound and unbound particles. By assuming that bound-states are meta-stable and close-to-equilibrium (implying that their number density Y B is almost constant, dY B /dx ≈ 0) one can describe the convoluted system of Boltzmann equations in an effective manner [27]. This assumption is reasonable as long as chemical equilibrium is assured between free particles (a premise we already made) and that the processes involving the bound-states (formation, ionization, rapid decays, leveltransitions, etc.) are fast enough to exceed the Hubble expansion, which is indeed the case for temperatures larger than their binding energies (see also Ref. [63] for a recent more detailed argumentation).
In a very similar fashion as for the bare co-annihilation scenario, we can now consider a single Boltzmann equation for the total number density of the dark sector particles (cf. Eq. (2.3)). The effects introduced by BSF can in fact be reabsorbed into an effective thermallyaveraged BSF cross-section affecting the evolution of the X number density and, thus, the total one, as shown in [27]. This effective contribution is given by where angular parentheses indicate thermal averaging.
We notice that, at large temperatures the ionization processes are much more efficient than bound state decays, hence the effective contribution of bound states in the dark sector evolution is negligible at these temperatures so that the relic density calculation is actually independent on the BSF cross-section. As the universe cools down, the ionization rate becomes exponentially suppressed and, eventually, decay processes will dominate, efficiently depleting the dark sector scalars, so that the effect of BSF on the Boltzmann equation will be potentially relevant even at temperatures close to the bound state binding energy (T E B ).
Since in our model we consider bound states formed by three (six) types of colored scalar pairs in the u R and d R (q L ) models, the BSF contribution will be given by the sum of three (six) terms like Eq. (3.24). In this spirit, we can write the total annihilation cross-section for the X i particle of flavour i by adding these terms to each σ X i X † i v rel . Summing over all the flavours, we obtain the following total effective annihilation cross-section for the X color triplet scalars: [1] . (3.25) This quantity supersedes the naive XX † annihilation cross-section in the first term in Eq. (2.11). Again, the importance of the new term becomes relevant at later times, where BSF can efficiently deplete the relic density way beyond the typical scales of thermal freeze-out via late decays of the dark sector bound states (see, e.g., Fig. 6 in [34]).
At this point, we would like to stress a non-trivial consideration. Given the fact that the color triplet scalars forming the bound states are unstable against Yukawa-mediated decays into DM and a quark, one could potentially encounter the situation in which the lifetime of a constituent particle is comparable or even shorter than the bound state lifetime itself, for example for large values of g DM . In this situation, the bound state could decay to B(XX † ) → X † +(χ+q) and this eventuality needs to be included when writing the Boltzmann equation. Although this additional decay channel might appear to play only the role of a bound state destroyer, thereby reducing the effect of BSF in depleting the DM abundance, by exploiting the principle of detailed balance in equilibrium, also the opposite reaction X † + χ + q → B(XX † ) must be considered. This sort of bound-state inverse decay provides an additional competing formation channel: its contribution can be recast into an effective BSF cross-section, with a form equivalent to Eq.(3.24). We refer the reader to Appendix B, for a treatment of the problem. Effectively, by assuming that the constituent-decay of the bound state is regulated by g DM , we show in Eq.(B.11) that the effective cross-section must be a monotonically-increasing function of the coupling g DM . In this sense, by neglecting the constituents decay, one is effectively considering the most conservative estimate of the effect of bound states on the DM abundance evolution. We leave a more detailed analysis of this effect for future work.
Notice that, in principle, bound states between colored scalars of same and different flavor i and j could exist, since their formation does not depend on the flavour structure. This is the case, for example, for a B(X i , X j ) particle-particle bound state, which however falls in the class of same-representation bound-states that have zero BSF amplitude in the degenerate limit (cf. Eqs. (A.8)-(A.10)). However, due to the flavor-diagonal structure of the Lagrangian considered in Eq. (2.2), different-flavour bound states, including particle-antiparticle bound state, cannot decay into gluons, but only into a quark-antiquark pair. Nevertheless, as Tab. 1 illustrates, this decay is less efficient than the decay into gluon-pair. First, for most of the parameter space the annihilation cross-section into quarks is velocity-suppressed. Second, even when considering the top-quark mass, the velocity-independent part of the perturbative cross-section XX † → qq that would determine the decay rate Γ B→qq becomes larger than the velocity-independent part of the XX † → gg only when g DM is large 14 . Therefore, the decay rate into qq can be safely neglected. The remaining decay channels that do not depend on the DM coupling (e.g., into γγ, Zγ or ZZ), one can see that they account for less than 1% of the total decay width of the bound state.
for the parameter space of interest, we do not consider the contribution of B(X † i , X j ) particle-antiparticle to the effective annihilation cross-section, given in Eq. (3.24).

Impact of the long-range effects on the t-channel model parameters
In the following, we investigate the impact of the Sommerfeld effect and bound state formation on the parameter space of the simplified t-channel models. In order to solve the Boltzmann equation, we employ the publicly available program micrOMEGAs 5.2.7 [52] and modify it to account for non-perturbative long-range effects on the thermally-averaged cross-section. Our code extends [35] by accounting for the emission of a gluon from the mediator (cf. Fig. 3) that was previously neglected in the literature and using the strong coupling constant according to the scale involved. Moreover, our implementation encompasses the depletion of the DM relic density by bound states that form and decay way beyond the typical time scales of thermal freeze-out of x 30, a feature not included in [35]. One can indeed check (see, e.g., Fig.6 in [34]) that the maximum of the bound state contribution to the effective annihilation cross-section can lie at temperatures of x 1000, allowing for DM relic density depletion at x 30. Finally, by making use of the assumptions illustrated in Sec. 3, our implementation fully exploits the ability of micrOMEGAs to automate the calculation of cross-sections by using the CalcHEP package [64], without the need to externally calculate them 15 .
In Fig. 4, we show the impact of non-perturbative effects on the parameter space of the u R , d R and q L models described in 2. For a fixed value of the Yukawa coupling g DM , each band represents the parameter range that is in agreement with the observed relic density of DM within 5σ acceptance Ω DM h 2 = 0.120 ± 0.005 [1]. The two different colors refer to two 14 For the parameter space here considered, one can check that, within the co-annihilating region, ∆ 0.2 mDM, in order to ensure ΓB→qq ΓB→gg one would need gDM > 2. For these values, we would always find under-abundant DM (e.g., cf. Fig. 5). 15 The code will be made available in an upcoming work [45].  different benchmark values for g DM , 10 −2 (orange) and 1 (blue). The different style of the lines identify the type of interactions considered: perturbative cross-section only (dotted), with Sommerfeld effect (dashed) and with Sommerfeld+BSF (solid). We notice, as expected, that non-perturbative effects are more pronounced for smaller values of g DM and smaller mass splittings ∆m. Table 1 shows that for a smaller ∆m/m DM the suppression of processes involving two colored scalars (cf. Figs. 1b-1g + BSF) is lifted such that BSF and processes that are affected by the SE become more relevant. Additionally, BSF is purely governed by the strong gauge coupling and thus less relevant if the competing processes are enhanced by a large g DM . 16 Therefore, BSF and the SE have an larger impact on the viable parameter space for a smaller g DM and for smaller ∆m/m DM . When fixing g DM , the total effective annihilation cross-section of DM is maximized for ∆m → 0, since the contributions from co-annihilations and colored annihilations are maximized. Thus, we find the maximal DM mass in agreement with the observed relic density for ∆m → 0. In the u R model, for small g DM , the maximal dark matter mass is shifted from roughly 1.0 TeV to 2.4 TeV. For larger couplings, g DM = 1, the maximal dark matter mass is expected at 2.9 TeV instead of 2.0 TeV.
For the d R model, the overall behavior is very similar and entails a shift of the upper bound on the dark matter mass from 1.2 TeV to 2.4 TeV for the smaller g DM and from 1.9 TeV to roughly 2.8 TeV for g DM = 1. Finally, in the q L case, the effect of BSF on the relic density is milder, as already anticipated at the beginning of this section, with a shift from 1.0 TeV to 2.0 TeV (g DM = 10 −2 ) and from 1.8 TeV to 2.5 TeV (g DM = 1.0), respectively. This effect is caused by the effective co-annihilation cross-section (see, e.g., Eq. (2.12) for the ∆m = 0 case). With more co-annihilating species involved, the effective contribution of each channel is smaller as illustrated in section 2 (cf., e.g., Eq. (2.12)). In the following, for conciseness, we will restrict ourselves to discuss the u R model, since the basic features of the d R and q L cases are very similar.
In Fig. 5, we display the values of g DM in the u R model that are able to account for the maximally viable observed DM density parameter Ω DM h 2 = 0.125 (at +5σ) in the (m DM , ∆m) plane, when BSF and Sommerfeld corrections are taken into account. In this sense, the displayed couplings serve as a lower bound on the Yukawa coupling g DM , since a smaller coupling would lead to an overproduction of DM.
We obtain a more stringent lower bound on g DM for an increasing DM mass, since a larger DM mass implies a smaller relic abundance, and therefore larger effective annihilation cross-section, in order to keep the density parameter at a constant value.
This behavior does not apply for points in the parameter space where the relative masssplitting is larger than roughly 20 %. In fact, for larger relative mass splittings, annihilation processes involving one or more colored scalars (cf. Figs. 1b-1i) are not efficient due to their 16 In fact, the strong coupling constant gs = √ 4παs lies in between 1.0 and 1.6 in the region of interest, depending on the average momentum transfer considered. Thus, for gDM = 1 interactions regulated by gs are comparable in strength to interactions regulated by gDM. Here, we take into account BSF+Sommerfeld processes, which become more relevant the closer we are to the co-annihilating region (small masssplitting). The dashed white lines are in correspondence of the region where 10 −7 g DM 10 −2 . Below this interval, it is not possible to ensure chemical equilibrium between the unbound particles in the dark sector and the co-annihilation freeze-out assumption breaks down.
Boltzmann suppression. This change of mechanism can be recognized by the sudden change in the bands in the small mass regime around the region corresponding to ∆m 0.2 m DM . Furthermore, decreasing the mass splitting implies a smaller g DM , as not only DM annihilations are less suppressed by the (smaller) mediator mass but, additionally, the exponential suppression of the colored annihilations is lifted; these latter are, in turn, enhanced by the Sommerfeld correction and BSF. Eventually, we encounter a region of parameter space where annihilations mediated by the strong gauge coupling g s are (almost) efficient enough to deplete the DM density sufficiently on their own. This region is enclosed by the white dashed lines and the lower bound on g DM lies in the interval 10 −7 g DM 10 −2 .
For even smaller mass splittings (in the gray region below the lower white dashed line), we are not able to find a lower bound on g DM anymore, because the corresponding coupling does not suffice to establish chemical equilibrium within the dark sector, so that the freezeout production of DM does not apply. In this region, DM freeze-out leads to a relic density (much) smaller than observed in the measurement of the cosmic microwave background. 17 We estimate the smallest couplingg DM ensuring chemical equilibrium in the dark sector by demanding that the interaction rate Γ X of the (inverse) decays X ↔ χq fulfills for T m χ /30. This condition ensures that chemical equilibrium between X and χ is established at the time of thermal freeze-out 18 . When neglecting the quark mass, assuming a small mass splitting ∆m and only considering two-body (inverse) decays of X, we obtain from Eq. (4.1)g DM m DM GeV 1 · 10 −9 + 6.8 · 10 −11 m DM ∆m , (4.2) Note, that Eq. (4.2) only applies when m DM ∆m m q and thus does not describe the chemical equilibrium condition of the mediator X associated with the top and DM where conversions proceed at next-to-leading order. We conclude that DM production does not proceed via thermal freeze-out for couplings g DM g DM .
In order to demonstrate the error on the lower bound on g DM (cf. Fig. 5) made when not including the SE and BSF, we compare in Fig. 6, the corresponding lower bound on the coupling g SE+BSF DM obtained when including BSF and SE with the lower bound on g Pert DM obtained considering perturbative annihilations only. It is important to note that the correction can lead to higher and lower values such that a dedicated analysis is necessary in order to estimate its implication on the parameter space.
Blue regions in Fig. 6 imply that a smaller g DM is expected in these scenarios when including SE and BSF, while orange regions indicate the opposite. The SE can, depending on the sign of the potential, enhance or decrease the effective annihilation cross-section, while BSF, which effectively provides a new annihilation channel, can only increase the effective annihilation cross-section of DM. Furthermore, the SE provides a correction to any process featuring two colored particles in the initial state (cf. Figs. 1b-1g), even if they are exclusively mediated by g DM . Thus, the SE has an effect even in the case g DM g s . BSF, on the other hand, as a new annihilation channel purely mediated by g s , becomes less important as soon as g DM g s . With that in mind, and by comparing Figs. 6 and 5, we find that orange regions, where the coupling required to not overproduce DM gets enhanced, coincide with regions of the parameter space where g DM g s ≈ 1. As we can extract from Table 1, this region of parameter space is dominated by direct DM annihilations (cf. Fig. 1a) and colored annihilations involving a XX initial state (e.g., cf. 1g), which results in a repulsive potential. For parameter points with ∆m 0.2 m DM , colored annihilations are strongly Boltzmann suppressed and non-perturbative effects are irrelevant. Blue regions, where the coupling required is reduced when considering the SE and BSF, are a result of strong bound state effects due to a small relative mass splitting ∆m/m DM and a smaller g DM g s . Furthermore, for g DM g s , colored annihilations involving a XX † −pair mediated by g s (cf. Fig. 1b-1f) dominate the effective DM annihilation cross-section, resulting in a attractive potential. While the effects induced by SE for a repulsive potential are on the percent level, the corrections to the coupling for regions with both SE for an attractive potential and efficient BSF are significantly more sizable. Due to this non-trivial behaviour throughout the parameter space, for a definite exclusion or discovery of the model, BSF and SE have to be considered. We want to emphasize that a constant correction factor ("K-factor") is not sufficient.

Limits on g DM from Direct Detection
To calculate direct detection constraints on the parameter space, we follow Ref. [37] 19 . Direct detection (DD) constraints arise from the non-observation of DM-nuclei scattering on earth.
The constraints on the DM-nucleon cross-section come from spin-independent (SI) and spindependent (SD) interactions. We use current spin-independent limits from Xenon-1T [4] and spin-dependent limits from the PICO-60 experiment [5]. Future projections are considered for the planned DARWIN experiment [67]. In our model, SD DM-nucleon scattering is mediated at tree-level by the s-channel exchange of a colored mediator X and the SD DM-nucleon crosssection increases with g 4 DM . For SI scattering however, due to the Majorana nature of the DM candidate, the velocity unsuppressed tree-level contribution is absent. Thus, SI DM-nucleon scattering is induced at the one-loop level, where it receives its dominant contribution from the diagrams shown in Fig. 7. Just as in SD scattering, the parametric dependence to the Yukawa coupling for the SI DM-nucleon cross-section also scales as g 4 DM . To compute the spin-independent DM-nucleon scattering cross-section in this simplified model, we perform a complete one-loop matching of the relevant Wilson coefficients. Hereby, we consider all possible diagrams and interference effects. We also perform a renormalization group evolution (RGE) from the scales of the mediator mass to the low-energy scale ( 1 GeV), relevant for DM scattering with the heavy nucleon. A detailed account of this is provided in [37]. Including the RGE evolution leads to an enhancement of roughly a factor of two at the amplitude level, hence roughly a factor of four at the cross-section level. Therefore, since the experimental limits on SI cross-sections are about six orders of magnitude stronger than for the SD ones, and since the SI cross-section, although one-loop-suppressed, is enhanced by the RGE-evolution, both SI and SD constraints are relevant in our study.
In Section 7, we will illustrate the impact of the DD constraints on the parameter space of the u R model here examined. Given the dependence of both SI and SD cross-sections on g 4 DM , from DD experiments we generically obtain upper limits on g DM for each data point (m DM , ∆m). This upper limit combined with the lower limit on g DM from preventing an overabundance of DM, will set the allowed and excluded areas of the parameter space, as we will discuss in detail in Ref. 7.1.

Collider Constraints
Since the mediator and dark matter have masses around the TeV-scale, the class of models considered are also subject to constraints from collider experiments. There are three possible searches that can be used to place constraints on the parameter space: 1. Prompt searches: here the mediators are produced and decay promptly to dark matter and jets.
2. Long-lived particle (LLP) searches: here the mediators that are produced in at the LHC do not decay promptly and may decay inside or outside the detector and require a dedicated search different from the prompt decay searches.
3. Bound state formation at LHC: here, mediators that are pair produced can form bound states and subsequently decay to SM gauge bosons.
We will show that prompt searches, LLP searches and bound state searches are relevant and complement each other in constraining the parameter space of the model. In the following subsections, we describe each of these LHC searches and our procedure for recasting them for the t-channel model.

Prompt Searches
Representative Feynman diagrams for the most relevant processes are shown in Fig. 8 and include: 1. pair-production of mediators followed by decay to a quark and dark matter.
2. Associated production of dark matter with a mediator.
3. Pair-production of dark matter with an Initial State Radiated (ISR) jet.
The production cross-section of the first of these processes contribution purely mediated by the QCD strong coupling, while the latter processes are controlled by the product of dark X g X † g X q g q X χ q g χ χ q X Figure 8: Feynman diagrams for processes at the LHC. The lines follow the same notation as in Fig. 7. In particular, from left to right, we draw the pair-production of mediators, the associated production of DM with a mediator, and the DM pair production with ISR jet.
matter coupling g DM and the strong coupling constant. Finally, processes scaling with g 4 DM such as pair production of DM and qq → XX are also included in the analysis. However, for the region of interest in this paper 20 , the phase space suppression ensures that these cross sections are also subdominant. The cross-sections for these processes are calculated at the Next to Leading Order (NLO) using MG5 aMC@NLO [68]. We use the CT14NLO [69] PDF set and set the renormalization and the factorization scale at µ r = µ f = H T 2 , where H T is the scalar sum of transverse momenta of all final state particles 21 . The details of the cross-section calculation is provided in [37]. As detailed in [37], the K-factors for these processes are fairly flat and range between 1.5-1.6.
Having described the Monte-Carlo procedure for the generation of events for the collider analysis, we now look at the relevant searches. Given the processes relevant for LHC referred to in Fig. 8 and enumerated above, prompt searches at the LHC constrain a bulk of the parameter space, as will be shown later. Both mono and multi-jet + E T searches at the LHC are relevant, especially for the region of parameter space where the mass difference (≥ 5 GeV) and couplings are sizeable. The mono-jet + E T searches primarily constrain the compressed spectra, and relies on an emitted hard ISR jet, while the multi-jet + E T search constrain moderate to high ∆m region, due to increased jet activity in the final state.
We will first start by describing the constraints arising from mono and multi-jet + E T searches at the LHC. They are obtained by recasting the most updated searches at the LHC for these processes, and reinterpreting them in our model. We use the MadAnalysis5 Public Analysis Database (Ma5-PAD) framework [71] to recast and reinterpret ATLAS searches. While we only use the latest ATLAS searches to obtain our constraints, we point out that CMS searches have similar coverage, and the constraints only differ slightly depending on the process. We briefly describe the validation process for each of the searches. The impact of these searches on the relevant parameter space in question for this work is discussed in Sec. 7.

Mono-Jet + E T Searches
We briefly summarize the recasting of ATLAS search for mono-jet + E T final state, described in [72]. In addition to trigger and isolation criteria of objects like leptons and jets, described in [72], the search requires the following pre-selection criteria: • A leading jet with p T > 150 GeV within a pseudorapidity gap of |η| < 2.4, and up to additional three jets with p T > 30 GeV and |η| < 2.8.
• A minimum missing transverse energy of 200 GeV.
• A requirement of the azimuthal angle between jets and missing energy of ∆φ(jet, E T ) > 0.4 for E T > 200 GeV and ∆φ(jet, E T ) > 0.6 for ( E T > 250 GeV), respectively.
• Isolated leptons are vetoed with a pre-selection criteria of p T > 20 GeV for electrons and p T > 20 GeV.
In addition to the pre-selection criteria, the search is split into 12 exclusive and 12 inclusive search regions depending on the missing energy criteria. We recast the analysis within the Ma5-PAD framework [71]. The recasted code, as well as the details of the recast process will be made public at [73]. A brief description of the recast process following the experimental search is described here. Jets are reconstructed using FASTJET [74], using the anti − k T [75] algorithm with a jet radius of 0.4, and with detector effects simulated by Delphes3 [76]. All parameters are tuned to match specifications provided by ATLAS. The analysis was validated by reproducing the cutflows for the benchmarks provided in the ATLAS documentation [72]. The cutflows were reproduced to within 10 % of the official documented results, and therefore we consider this analysis to be validated. For the re-interpretation purpose, signal events corresponding to the processes listed above were generated at tree level in MadGraph5 [68] with up to two additional jets. These signal events were subsequently showered and hadronized using PYTHIA8 [77], with matrix element and parton shower (ME-PS) [78] merging, with a merging parameter set to m X /4. The event rates after all cuts are normalized to the NLO cross-sections and an integrated luminosity of 139 f b −1 . To extract the limits on the models, based on the normalized number of events obtained for each signal region, we apply the log-likelihood method to obtain an exclusion for each point in the m χ − ∆m plane. An in-built confidence level calculator within Ma5-PAD was used to compute the associated upper limit at the 95 % confidence level (CL) on the signal cross-section according to the CLs prescription [79]. Even though the analysis contains a large number of signal regions, the bulk of the exclusion limit is determined predominantly from signal regions that have large signal rates, low background rates and small uncertainties.

Multi-Jet + E T Searches
The class of multi-jet + E T searches targets BSM scenarios with at least two jets and a significant amount of missing transverse energy. We use the latest ATLAS search for this channel at an integrated luminosity of 139 fb −1 for re-interpretation purposes [70]. The search targets SUSY particles, in particular gluinos and electroweak gauginos within a simplified model scenario that undergo cascade decays to yield a large multiplicity of jets and E T in the final state. We re-interpret this search to set constraints on the models described in this paper.
The baseline criteria for this search is, • At least two isolated jets with p T > 20 GeV and |η| < 2.8.
• A leading jet with p T > 200 GeV and a subleading jet with p T > 50 GeV.
• A criteria of the azimuthal angle between the leading two jets and missing energy ∆φ(jet, E T ) > 0.2 • An effective mass (sum on the scalar p T of all jets and E T ) m eff > 800 GeV.
In addition the search is divided into a large number of signal regions depending on the final-state targets and mass gaps. For search regions without final states involving b-jets, which is what we require for our re-interpretation purposes, the signal regions are designed to optimize the signal-to-background ratio in bins of p T , the number of jets and E T . The details of the full search strategy can be found in [70]. To recast this analysis, we follow the experimental paper, with the implementation done within the Ma5-PAD framework. Pre-selection criteria, detector specifications and signal regions were implemented within the fast simulation platform. For benchmark points documented in [70] for the recast were generated using MadGraph5 [68] with at least two additional jets. Parton showering and ME-PS matching was performed with PYTHIA8. As was the case with the mono-jet + E T search, jets were constructed with FASTJET with a radius of 0.4 using the anti − K T algorithm. The signal yields, and the subsequent cutflows were compared with the official ATLAS documentation. We find that the recasted search agrees with the official cutflows within 10 %, and hence we consider the analysis to be validated.
To obtain the projections for high-luminosity LHC (HL-LHC) [80], we provide a naive luminosity-rescaled 95 % confidence level exclusion limit. In order to do this, we simply rescale the signal and background yields to HL-LHC in order to estimate the exclusion limit. We note that these are rough estimates, while more accurate limits depend on a proper background estimate from data-driven methods, dealing with high pile-ups as a consequence of a larger bunch crossing at high luminosity.

Bound State Production at LHC
Colored mediators can be pair produced at the LHC and subsequently form bound states B(XX † ). When the coupling of the mediator to dark matter (g DM ) is small, these bound states predominantly decay to pairs of gauge bosons. For the specific scenario we consider here, the decay to pairs of gluons is dominant, with a branching ratio ∼ 99%, followed by decay to pairs of photons, Zγ and then pairs of Z-bosons. 22 Here, we consider the production of the bound state B(XX † ) and its subsequent decay to pairs of photons at the LHC, since the diphoton decay channel provides the strongest constraints for the u R model. Importantly, these constraints do not depend on the mass of dark matter and, for small values of g DM , are not sensitive to the mediator and dark matter coupling.
We therefore use results from high-mass ( 100 GeV) diphoton searches performed by the ATLAS experiment [81] to place constraints on the mass of the mediator. The analysis carried out in Ref. [81] provides limits on the fiducial cross-section (times branching ratio) at 95% C.L., which we correct by an almost flat efficiency factor, given by the detector acceptance times the reconstruction efficiency, as described in [82], in order to obtain the corresponding total cross-section (times branching ratio). In order to calculate the theoretical cross-section, we follow the procedure described in Refs. [83,84], which calculates the boundstate production cross-section at NLO (QCD) and places constraints on the s-wave (J P C = 0 ++ ) stoponium production at LHC. In the following, we briefly outline the main key points of the analysis. The LO production cross-section for a stoponium-like bound state in the Narrow Width Approximation (NWA) is given by Here, Γ B(XX † → gg) is the decay width of the bound state to pairs of gluons, P gg is the gluon luminosity for proton-proton collisions at 13 TeV of center-of-mass energy, and m B 2m X , is the mass of the bound state 23 . We exploit the NLO calculation derived in [83] and evaluate the strong coupling α s at the scale µ = m B ; we then multiply this result by the branching ratio of bound state decays into diphoton resonances (around 0.3% − 0.5%). We refer the reader to [83,84], where relevant expressions for the decay width of B(XX † ) to gg, γγ, Zγ and ZZ can be found. Accounting for the three generations of stoponium-like bound states of the u R model, we show in Fig. 9 the present exclusion limits on the bound state mass by comparing the experimental bound on the production cross-section times the diphoton-resonance branching ratio (black line) to the theoretical one (red line) for the u R model.  [81]. The red line is our theoretical prediction for the u R model considered. We exclude masses whenever this line is above the experimental one.
the exclusion limits on the mediator mass m X m B /2: 100 GeV m X 290 GeV . (6. 2) The lower limit of 100 GeV is an artifact of the experimental analyses and lower masses can be probed by looking at data from other lower energy experiments.
In order to provide an estimate of the projected exclusion bounds from the HL-LHC, we perform a rescaling of the integrated luminosity with respect to the data used from ATLAS previously 24 . The upper bound on the mass of the mediator increases significantly for HL-LHC and we find that m X 650 GeV .
We would like to stress once again that these limits are independent of the precise value of g DM , the Yukawa coupling connecting DM with the mediator and SM particles. 24 We use the asymptotic expressions to determine the significance (Z = signal/ √ background). The number of background events are determined from Fig. 2 of [81] and are rescaled by the luminosity for HL-LHC. We cross-checked that this method is consistent with the experimental bounds obtained from the full experimental likelihood analysis, as has also been noted in [82]. Additionally, we checked that HL-LHC will have a large number of background events so that the asymptotic form of the significance is a good approximation.

Long-Lived Particle Searches
A region of the parameter space of our interest with small mass gaps results in decay widths that can only be probed by LLP searches. For this work, we consider searches for heavy stable charged particles (HSCP).
Pair-production of mediators result in charged particles, that, for certain regions of parameter space, can have small enough decay widths such that they decay outside the detector. Since the mediators we study in this model are color triplets, they will form R-Hadrons (neutral or charged). The exact dynamics of the formation of charged hadrons is governed by QCD. For the purposes of this work, we will use the cloud model of hadronization referred in [85]. The heavy charged hadron with a velocity β = v/c < 1 travels through the detector and decays after crossing the tracker or the muon chamber, leaving an ionizing track with an energy loss larger than SM particles. For R-hadrons that decay outside the detector, the time-of-flight (TOF) measurement can distinguish BSM particles from the corresponding SM particles.
ATLAS and CMS have performed HSCP searches at 8 and 13 TeV center-of-mass energies at the LHC, with results presented within supersymmetric models with long-lived gluinos and squarks. We follow the method employed in [86] to recast and re-interpret CMS searches for HSCP at 8 and 13 TeV LHC [85,87]. The CMS search constrained long lived stops, staus and heavy vector-like leptons utilizing the tracker only and the tracker+TOF analysis. Note that the tracker+TOF analysis is relevant only for large values of cτ (≥ 10 m). We directly translate the cross-section limits on stops and staus obtained by CMS to constraints on our model. The regime of maximal sensitivity depends on the lifetime of the parent particles. For cτ 1 cm, a significant fraction of LLPs decay within the tracker, resulting in a suppression of the HSCP signal. Therefore, we multiply the pair production cross-section by the fraction of particles decaying either outside the tracker or outside the CMS detector. This fraction is dependent on trigger and selection criteria. Following previous works [86], we obtain the effective cross-section by multiplying the efficiency fraction f LLP (τ ), detailed in [88] as, where L is the detector size, which corresponds to L = 3(11) m 25 for the tracker only (tracker+TOF) analysis. The effective cross-section is then directly compared to the CMS cross-section upper limits for direct production of stops. The cross-sections at leading order is calculated using MG5 aMC@NLO [68]. We finally provide a tentative projection for HL-LHC following previous work in [86]. As noted there, high-luminosity LHC projections for LLP searches get complicated by the fact that backgrounds are primarily instrumental and therefore cannot be simulated using Monte-Carlo generators. Moreover, HL-LHC requires new triggers as well as a better understanding of pile-up rates. Nevertheless, following [86], we perform a simple re-scaling of expected signal and observed background events. To this end, we normalize the signal events obtained after all cuts to the production cross-section and an integrated luminosity of 3000 fb −1 . For the background, we take the official number of background events from the CMS analysis and normalize it for HL-LHC.

Results and discussion
In the following, we demonstrate how current and future experimental measurements constrain the parameter space of the u R t-channel simplified model and how the inclusion of Sommerfeld effect (SE) and bound-state formation (BSF) affect these exclusion limits. We show the resulting exclusion limits in Fig. 10, obtained by considering only perturbative annihilations (a), with the inclusion of the SE (b), and by considering the SE as well as BSF (c), respectively. In order to arrive at these limits, we first determine the smallest coupling g DM,Ω that does not overproduce DM for a given (m DM , ∆m) point (cf. Fig. 5), and then compare it to the upper limits on g DM from direct detection (DD) and collider experiments g DM,exp . If g DM,Ω > g DM,exp , we regard the data point as excluded by the experiments considered.
While the precise value of g DM,Ω is determined by a scan with a modified version of micrOMEGAs, as described in section 4, the scaling behavior of g DM,Ω can be understood by means of a simple estimate for the relic density, Ω DM ∼ σv −1 . For this estimate, it is useful to distinguish between three different regimes. Firstly, the region of parameter space where DM annihilations are mainly set by colored annihilations, which arises for DM masses below and around a TeV while ∆m m DM . Thus, g DM,Ω can be chosen almost arbitrarily small, as long as chemical equilibrium in the dark sector is assured, such that the lower bound on g DM,Ω approaches zero. Secondly, the region of parameter space where colored annihilations efficiently contribute to σv but annihilations mediated by g DM are necessary to sufficiently deplete the relic density and provide the dominant contribution to σv . This happens for larger DM masses m DM O (TeV) and ∆m m DM . Lastly, the region where colored annihilations are suppressed by a large mass splitting ∆m m DM and the annihilation cross-section is purely set by direct DM annihilations. This estimate results in In the following, we comment in detail on the implications of the different constraints.

Interplay of current experimental constraints
Spin-Independent Direct Detection (SI-DD). Following the discussion in Section 5, we show in Fig. 10, a green area corresponding to the resulting excluded regions from SI-DD searches in the (m DM , ∆m)-plane. They lead to strong exclusion limits for small mass splittings (∆m 100 GeV) and large DM masses (m DM 1 TeV). It is noteworthy,  For each point in the plane, we determine the smallest value of g DM such that DM is not overproduced. If this lower bound on g DM contradicts the limits from spin-independent DD (spin-dependent DD, prompt-collider searches, perturbative unitarity, stoponium searches), it is colored in green (magenta, blue, black, cyan). We show the experimental limits obtained considering perturbative annihilations, the SE and the SE+BSF. The gray dashed line divides the parameter space into two regions. Above, the observed relic density can be generated via thermal freeze-out. Below, thermal freeze-out under-produces DM and accounting for the complete DM relic density requires production via conversion-driven freeze-out or freeze-in.
that SI-DD is unable to constrain mass splittings ∆m 50 GeV for smaller DM masses, m DM 1 TeV. The size of this region strongly depends on the inclusion of non-perturbative effects, as we discuss in the following section (cf. Section 7.2). Additionally, we observe an excluded region around mass splittings equal to the top quark mass, which is a result of a top-quark resonance in the DM-nucleon scattering cross-section calculated at one-loop order [37]. The scaling behavior of the SI-DD constraints can be understood directly from the form of the DM-nucleon cross-section [37] σ SI ∝ where m q is the mass of the quark propagating in the loop mediating the interaction of DM with gluons and m N is the mass of the nucleon. Additionally, the lower bound on g DM is fixed by the requirement of not overproducing the relic density. Using Ω DM ∼ σv −1 we arrive at Eq. (7.1). Inserting Eq. (7.1) into Eq. (7.2), we obtain This result supports the observations made previously. Spin-independent DD is able to constrain areas of the parameter space involving small mass splittings that are not dominated by annihilations mediated by the strong gauge coupling, whereas regions of large mass splitting remain unconstrained due to the suppression by six powers of the relative mass splitting.
Spin-Dependent Direct Detection (SD-DD). Following the discussion in Section 5, we depict the resulting SD-DD excluded regions with magenta tones in the (m DM , ∆m)-plane in Fig. 10. Spin-dependent DD is more constraining than spin-independent DD for larger mass splittings ∆m 200 GeV, while the exclusion limits are less stringent for mass splittings of ∆m 100GeV and, as for SI-DD, we find an unconstrained region of parameter space in the case of m DM 1 TeV and ∆m 50 GeV. Again, this fact can be understood from the scaling of the SD-DD cross-section: where in the second step we inserted the estimates for the lower bound on g DM given in Eq. (7.1). Comparing the results for the SD and SI DM-nucleon cross-section reveals that for ∆m m DM the SD contribution is only enhanced by two powers of the relative mass splitting and thus less constraining than the SI contribution, which is enhanced by four powers of ∆m/m DM . On the other hand, for large mass splittings, the SD contribution does not depend on ∆m and thus is able to constrain these areas of the parameter space below a certain DM mass threshold of m DM 200 GeV.
Partial-wave unitarity bounds. In order to ensure the perturbativity and unitarity of the model, we require a conservative limit of g DM < √ 4π. We show with a black shaded area in Figs. 10-12 the range of the parameter space when this condition is violated. The limit on g DM from partial wave unitarity for this simplified t-channel DM model was given in [89], relying on methods developed in [90].

Prompt Collider Constraints.
Following the discussion in Section 6, we show with blue colours the resulting excluded regions in the (m DM , ∆m)−plane in Fig. 10. They provide the most stringent constraints for m DM 1.2 TeV ∧ ∆m 100 GeV. We do not consider exclusion limits from prompt collider searches if g DM > √ 4π, since it violates the limit obtained from perturbative unitarity and additionally the assumption of a narrow width for the colored scalar X breaks down. For low mass gaps, the mono-jet constraints provide the strongest constraints as it relies on an emitted hard ISR jet recoiling against the missing energy. For large mass gaps, the jet multiplicity and the jet transverse momentum p T is significantly larger, and therefore the primary constraints originate from the multi-jet + E T . Overall, the difference between the two searches is not significant as the multi-jet + E T , this being an inclusive search, is generally competitive in most of the parameter space. The primary reason for this being the hard radiated jets from ISR that passes the somewhat strict pre-selection criteria for multi-jet searches. In summary, both multi-jet and mono-jet + E T searches yield similar limits in the low mass gap region, while the multi-jet + E T search yields stronger limits in the large mass gap region. Generally, the prompt searches have two effects that impact its shape in the (m X , ∆m)−plane. The first, in regions with small mass gap, the cross-section is enough such that there is significant sensitivity to existing searches. The second, regions of large mass gap, there is a combined effect of rising sensitivity and falling cross-section resulting in a cliff.
In addition, we provide limits from bound state searches at the LHC as described in section 6.1.3. We find that ATLAS excludes masses of the scalar mediators of m X = m DM + ∆ 290 GeV. The corresponding excluded data points are colored in cyan. In order to exclude a data point by this search, we further demand that the binding energy of the bound states exceeds the decay rate of the constituent, otherwise, bound state production at the colliders is suppressed [91]. We comment specifically on the potential of this type of search in Section 7.5.

Impact of Sommerfeld effect and bound state formation on the combined exclusion limits
Considering the experimental constraints from direct detection and collider searches as well as the unitarity bound, the remaining unconstrained parameter space can be divided into two distinct areas, requiring very different values of the coupling g DM (cf. Fig. 10). For large m DM and ∆m, it is required to have g DM O (1) in order to avoid the overproduction of the DM relic abundance; on the contrary, for smaller masses, the couplings lie in the range 0 ≤ g DM 0.1. Due to the considerable coupling strength g DM required, the first scenario is only mildly affected by SE and BSF.
As stressed in Section 4, BSF is negligible for larger DM masses and mass splittings. In this region, however, in the case with g DM > g s , the particle-particle colored annihilation cross-sections are Sommerfeld-affected by the presence of a repulsive effective potential 26 . Hence, a larger g DM is necessary in order not to overproduce dark matter such that the limits arising from direct detection, collider experiments and unitarity are stronger. The size of the effect is at the few % level and displays as a slight shift in, for instance, in the direct exclusion limits for large DM mass of 20 TeV in Fig. 10.
Remarkably, the parameter region featuring m DM 2.5TeV and ∆m 80 GeV is strongly affected by SE and BSF. To highlight this effect, we refer to the exclusion limits in Fig. 11. In this region, both the SE and BSF significantly enhance the annihilation crosssection of DM, thereby, allowing for smaller values of g DM . Consequently, the experimental limits are relaxed. For instance, we notice that the regions constrained by direct detection experiments are significantly reduced for almost all the mass splittings displayed in Fig. 11. Most notably, when including the SE and BSF, we find that parts of the parameter space for 500 GeV m DM 1.5 TeV, previously thought to be excluded, can indeed produce the observed relic density and not being ruled out by current experiments.
We observe that the largest viable DM mass in this region is shifted from m DM 1 TeV to m DM 2.4 TeV, while the maximally viable mass splitting increases from ∆m 30 GeV 26 The effect of Sommerfeld repulsion is only relevant when colored annihilations are most efficient, namely when ∆m 0.2mDM. As shown in Tab. 1, the only velocity non-suppressed process that depends on g 4 DM is XiXi → qiqi with same flavor i, whose color contribution results only in the sextet configuration, which is repulsive. As soon as gDM > gs this process dominates.  Figure 11: Exclusion limits from various experiments in the co-annihilating area. As before, the limits from spin-independent DD (spin-dependent DD, colliders, unitarity, stoponium searches) are colored in green (magenta, blue, black, cyan). In addition, we provide an estimate for the potential reach of LLP searches constraining the area of the parameter space where the correct relic density can be reproduced only via non-thermal mechanisms, such as conversion-driven freeze-out or freeze-in. Those regions are marked in orange.
to ∆m 50 GeV. In conclusion, considering the SE and BSF leads to a O (100%) correction on the exclusion limits in the strongly co-annihilating area. Finally, we would like to point out that only considering SE, although it provides a correction of O (10%) to the exclusion limits, is insufficient in the strongly co-annihilating area, as can be clearly observed from Fig. 11. Therefore, we want to emphasize that for a conclusive statement on the exclusion of such models, a full analysis including the Sommerfeld effect and bound state formation needs to be taken into account.
In the gray area shown in Figs. 10 and in the zoomed version Fig. 11, thermal freeze-out always leads to underabundant DM, which is not experimentally excluded but requires an additional mechanism to explain the observed relic abundance. For example, for couplings g DM small enough such that DM is not in chemical equilibrium with the colored scalars X i , DM production can proceed via alternative mechanisms such as conversion-driven freeze-out or freeze-in.
Finally, thermal freeze-out is able to account for the observed DM abundance in the remaining white regions. As clearly visible from Fig. 11, SE and BSF lead to a significant increase of the allowed parameter space, which would have been underestimated when considering only a pure perturbative tree-level calculation and hence these effects need to be included for a comprehensive study.
As a final remark, since the couplings required to generate the observed relic density below the gray dashed line are tiny, g DM 10 −6 , they can hardly be constrained by prompt collider signatures or direct detection experiments. However, for couplings of this size, the colored scalars X are potentially collider-stable and therefore they could appear through long-lived-particle (LLP) signatures. In what follows, we estimate the potential reach of LLP searches within these regions of the parameter space.

Potential of long-lived particle searches
In order to address the parameter space unable to account for 100% of the DM relic density via thermal freeze-out (and possibly tested by LLP searches), we estimate the sizeg DM below which the interaction allows for an out-of-equilibrium dark sector. To this end, we exploit the condition stated in Eq. (4.1), which ensures that the interaction rate of the (inverse) decays of the colored scalars is out-of-equilibrium at the typical timescales of thermal freeze-out. In this sense,g DM provides a rough estimate of the couplings involved in scenarios relevant for the conversion-driven freeze-out, while it certainly acts as an upper bound on couplings involved in freeze-in scenarios.
In contrast to the aforementioned constraints, the search for LLPs results in lower limits on the coupling g DM > g LLP DM , since a feebler coupling implies more stable colored scalars X i , which are therefore subject to tighter constraints from heavy stable charged particles (HSCP) searches, as described in section 6.2. We mark the areas wheng DM < g LLP DM in Fig. 11 with orange meshes: they indicate the potential of LLP searches to constrain DM production relying on an out-of-equilibrium dark sector 27 .
Overall, we find that HSCP searches could test DM masses up to roughly 1.2 TeV and decay lengths of O (10 cm). As a consequence, the complete area below the gray dashed-line is constrained when considering perturbative annihilations only. Still, a large part of it would remain unconstrained when SE and BSF are taken into account, since it would extend to DM masses beyond the current reach of HSCP searches. In conclusion, we stress that the inclusion of the SE and BSF is essential for determining the parameter space of interest for LLP searches.

Projected exclusion limits of future experiments
In Fig. 12, we illustrate the expected exclusion regions considering the projected sensitivities for the DARWIN experiment, as described in Section 5, and the HL-LHC, as discussed in Section 6. We find that, given the much larger amount of data expected for HL-LHC and considering the improvement of roughly two orders of magnitude on the excluded DM-nucleon cross-section with respect to XENON1T and PICO-60, the parameter space for freeze-out production of the observed DM relic density that was left unconstrained by the current experimental limits (cf. Fig. 10) will be tested almost completely (cf. Fig. 12). Notice that, while HL-LHC will enlarge the probed parameter-space, the most dramatic improvement will in particular arise from DD experiments.
In particular, assuming that the Majorana particles account for 100 % of the observed DM relic density, both SI and SD constraints will be able to almost entirely test the currentlyremaining parameter space of the u R model. Notice that, while Fig. 12 seems to imply that all the parameter space above the gray dashed line is excluded by direct detection, in fact a tiny white space is remaining as the green and gray dashed line do not overlap exactly. The remaining region corresponds to small couplings, 10 −6 g DM 10 −2 , where direct detection is evaded but DM production still proceeds via thermal freeze-out 28 . While the couplings involved are too small to allow for a signal in direct detection experiments or prompt searches, they are also too large to induce long-lived particle signatures. Remarkably, as visible in Figs. 11 and 12, bound state searches at the LHC manage to close this gap below a certain mass threshold. The latter is expected to increase significantly, from roughly 290 GeV at the LHC to roughly 650 GeV, at the HL-LHC.
For larger mass splittings, the experimental limits improve drastically, with SD constraints being especially powerful in this region. This is a result of the less pronounced 27 More precisely, we underestimate the exclusion regions for freeze-in, since g freeze-in DM <gDM. On the other hand, the limits only provide a rough estimate for scenarios of conversion-driven freeze-out, which is relevant at gDM ≈gDM. To arrive at a conclusive statement, a dedicated study of this region is required. For the conversion-driven freeze-out scenario an analysis of the dR-version of the model has been recently provided in [38]. The limits obtained therein are weaker since it is found that conversion-driven freeze-out can have sizable effects for couplings a factor 2 to 5 larger than obtained with our estimate in Eq. (4.1). Thus, in our work, we find decay lengths roughly a factor 10 larger than in [38], resulting in more stringent limits. 28 This region is narrow in ∆m (less or about a GeV), since the efficiency of the colored co-annihilations dominating the effective annihilation cross-section is regulated by exp 2∆m · m For each point in the plane, we determine the smallest value of g DM such that DM is not overproduced. If this lower bound on g DM contradicts the projected limits from spinindependent DD (spin-dependent DD, prompt-collider searches, perturbative unitarity, stoponium searches, LLP searches), it is colored in green (magenta, blue, black, cyan, orange). We show the experimental limits obtained considering perturbative annihilations, the SE and the SE+BSF. The gray dashed line divides the parameter space into two regions. Above, the observed relic density can be generated via thermal freeze-out. Below, thermal freeze-out under-produces DM and accounting for the complete DM relic density requires production via conversion-driven freeze-out or freeze-in.
suppression with the relative mass splitting for SD DM-nucleon scattering (cf. Eq. (7.4)) compared to the suppression of the SI scattering (cf. Eq. (7.3)).
If future experiments do not observe a signal, three areas of the parameter space will still survive. Firstly, a small region of freeze-out DM involving g DM ∼ O (1) couplings with DM and mediator masses in the multi-TeV regime. This scenario could possibly be probed by future collider experiments with larger center-of-mass energies such as the FCC [92]. Secondly, a narrow belt involving couplings 10 −6 g DM 10 −2 at the transition region to the non-thermal production area. Lastly, the region below the gray dashed line, where DM is either underproduced via freeze-out or created via conversion-driven freeze-out or freeze-in. Importantly, we find that the latter two regions are significantly affected when Sommerfeld effects and bound state effects in the early universe are taken into account. If DM production in this region is non-thermal, it can be tested by LLP searches. Again, we stress that LLP searches beyond the HL-LHC appear to be necessary to thoroughly probe this scenario.
To conclude, both in the absence or in the presence of a positive detection signal, it becomes evident that in order to accurately determine the remaining viable parameter space of simplified t-channel DM models, the inclusion of non-perturbative effects (SE and BSF) from long-range forces is mandatory when calculating the DM relic abundance.

Potential of Bound State Searches at the LHC
In order to highlight the potential reach of searches for dark sector bound states at the LHC (cf. sections 6.1.3 and 7.1), we show in Fig. 13 the parameter space exclusion limits in the (m DM , g DM )−plane for fixed values of the relative mass splitting δ = ∆m/m DM = 0.01 (upper panels) and δ = 0.05 (lower panels). The exclusion plots exploit the same color-scheme for the various experimental limits as in the (m DM , ∆m)-plane of Figs. 10-12. In contrast to the previous figures, it also features regions of parameter space where the DM abundance is overproduced, since the mass-splitting is fixed to the value ∆m = δ m DM , while the coupling is varied 29 . The portion of parameter space excluded due to overabundant DM is bounded by solid black lines and colored in gray. The correct DM abundance in the freeze-out regime coincides with the black solid boundary between the overabundant and freeze-out region 30 .
Direct detection and prompt collider searches are able to constrain the DM-SM coupling g DM from above, since the production rate at the LHC and the DM-nucleon coupling increase with g DM . Additionally, prompt collider searches require a minimal value for g DM to be efficient in order for the decay of the mediator to be prompt. On the contrary, long-lived particle searches can constrain the coupling g DM from below but cease to be efficient when g DM becomes larger, as the decaying particle would disappear too quickly to be tagged as an LLP. 29 The parameter space presented in Fig. 13 corresponds to a straight line with ∆m = δ mDM in the Figs. 10-12. In Fig. 13, however, we vary the coupling gDM, which in the Figs. 10-12 is fixed such that DM is not overproduced. 30 The correct relic density can potentially be obtained in the non-thermal region. However, we did not analyze this situation in detail.  The limits from spin-independent direct detection (spin-dependent direct detection, prompt collider searches, long-lived particle searches, bound state searches at colliders) are shown in green (magenta, blue, orange, cyan). Additionally the plane is split into three regions. Firstly a region for small couplings where DM production proceeds non-thermally, secondly a freeze-out region where DM is either underproduced or matches the observed relic density, which coincides with the boundary of this region to the area of parameter space where DM is overabundant.
The production cross-section for a X − X † bound state at a hadron collider is given by Eq. (6.1) and solely depends on the strong gauge coupling as long as the decay width of the constituent, regulated by g DM , is smaller than the binding energy of the bound state [91]. For larger values of the decay width, however, the production of bound states at colliders is suppressed since it is more difficult for a bound state to form before one of the potential constituents decays and thus bound-state searches are only effective up to a certain coupling, roughly when g DM g s . However, in contrast to prompt collider and long-lived particle searches, bound-state searches do not rely on a decay of the bound state mediated by g DM . Since the bound state confines a particle X and an antiparticle X † , it can decay via gaugemediated annihilations of the constituents, a scenario not viable for the decay of a single particle. As a consequence, bound state searches can constrain a large range of couplings g DM as long as it is small enough to allow for the efficient production of bound states at colliders. This can be clearly seen in Fig. 13, where the corresponding cyan region is the only search covering couplings in the range of 10 −6 g DM 10 −2 . Notably, the expected limit obtained via bound-state searches increases substantially at the HL-LHC, where it could potentially exclude DM production in the strongly-coannihilating region up to DM masses of 650 GeV.
We would like to stress that bound-state searches at colliders appear to be one of the most suitable search strategies to fully test the last remaining areas of the parameter space for colored co-annihilation scenarios featuring couplings of 10 −6 g DM 10 −2 .

Conclusions
We have studied the phenomenological implications of non-perturbative effects arising from long-range forces between color-charged dark sector particles, in particular the Sommerfeld effect (SE) and Bound State Formation (BSF) on the interpretation of experimental exclusion limits based on the example of a simplified t-channel dark matter (DM) model. Hereby, we analyzed a class of simplified t-channel dark matter models featuring a Majorana fermion χ, the DM candidate, coupled via a Yukawa interaction with SM quarks and a color triplet scalar X with the same quantum numbers of the quark it interacts with. In the strongly coannihilating regime, the DM abundance is mainly set by annihilations of the colored scalars, which are affected by the SE that can enhance or diminish the effective DM annihilation cross section depending on the corresponding potential. Additionally, bound states can form in the dark sector and their decays into SM particles provide another channel for depleting dark sector particles and hence DM. These effects are a natural consequence of the existence of color charge in such models and hence need to be taken into account.
We have shown that the size and the sign of the corrections to the DM annihilation cross section crucially depend on the coupling structure (g DM vs g s ) of the model (cf. Fig. 6). The SE affects any annihilation channel with two colored particles in the inital state, independently of the size of g DM . However, it enhances the annihilation cross section for g DM g s , while it tends to decrease the cross section otherwise, especially for g DM > g s . In contrast, BSF always increases the annihilation cross section since it acts as an additional annihilation channel. Being purely mediated by g s , BSF becomes irrelevant as soon as g DM > g s . Consequently, non-perturbative effects increase the DM annihilation cross section for g DM g s , while they can be irrelevant or even decrease the annihilation cross section otherwise. This clearly illustrates that a naive analysis of the non-perturbative corrections with simple flat factors is not an accurate solution, whereas a more detailed analysis of the processes and relevant model is crucial.
It is important to stress that when including SE and BSF, we predict larger dark matter masses than naively expected from tree-level calculations only. This is in particular interesting with respect to future indirect detection searches such as CTA. For instance, for a coupling of g DM = 10 −2 g s , we found that the consideration of SE and BSF leads to an increase of the expected DM mass m DM required in order not exceed the observed relic density. In the limit of a degenerate dark sector we predict as the maximal dark matter mass 2.4 TeV instead of 1.0 TeV (correction of around 140%) for DM couplings to right-handed quarks and 2.0 TeV instead of 1.0 TeV (correction of around 100%) for DM couplings to left-handed quarks, cf. Fig. 4.
Most importantly, we find that by considering SE and BSF, parameter space thought to be excluded by direct detection experiments and LHC searches remains still viable, cf. Fig. 10 and Fig. 11. This concerns especially the strongly coannihilating region, where nonperturbative effects have the largest impact. Based on tree-level calculations only, the viable parameter space was thought to span an area around (m DM , ∆m) (1 TeV, 30 GeV), however, when considering SE and BSF it opens up again to (m DM , ∆m) (2.4 TeV, 50 GeV). Additionally, we find that DM freeze-out can lead to the observed relic abundance while not yet being experimentally excluded for mass splittings ranging roughly from 50 GeV ∆m 3 TeV and DM masses ranging roughly from 200 GeV m DM 17 TeV for the u R model, apart for a region around the ∆m m t , where the one-loop-induced SI DM-nucleon scattering is resonant.
SE and BSF extend the parameter space which leads to an underabundant relic via the freeze-out mechanism, featuring very small couplings g DM ≈ 10 −7 − 10 −8 . We demonstrated that this parameter space can be probed by long-lived particle (LLP) searches, in particular heavy-stable-charged particle (HSCP) searches are able to probe ∆m 40 GeV and m DM 1.1 TeV, cf. Fig. 11.
In the strongly coannihilating region, thermal production of DM can generate the observed DM relic abundance within a narrow slice in the (∆m, m DM )-plane with couplings varying over four orders of magnitude (10 −6 g DM 10 −2 ). This coupling range is currently inaccessible to prompt collider searches, direct detection, as well as long-lived particle searches (LLPs). Remarkably, we found that searches for dark-sector bound-state resonances at colliders provide a complementary probe in this region of parameter space. This is due to the fact that prompt collider searches rely on large g DM necessary for a sizeable production cross section and the successive prompt decay of the colored mediator. This is in contrast to the production and decay of the bound state which proceeds via gauge couplings. Hence, bound state searches at the LHC are necessary to probe the complete parameter space relevant to colored coannihilation.
Finally, we investigated the reach of future DD and collider experiments, such as DAR-WIN and HL-LHC, and found that almost the entire remaining parameter space will be the tested in the near future (cf. Fig. 12). Only two small regions would still allow for thermal freeze-out as a viable solution: a corner around ∆m ∼ 2 − 4 TeV and m DM ∼ 3 − 6 TeV and an even narrower slice in the strongly co-annihilating regime, featuring tiny DM couplings. While the first one could be finally probed by colliders with higher center-of-mass energies (e.g., the FCC), the latter is potentially (and possibly uniquely) testable by dedicated bound state searches at colliders. This will be complemented by future HSCP searches at the HL-LHC, improving the limits in the underabundant/non-thermal region, c.f. Figs. 12-13).
In conclusion, in order to arrive at a final statement on the feasibility of the class of simplified t-channel models here analyzed to describe the observed DM relic density in the universe, non-perturbative effects arising from long-range forces must be taken into account. These are naturally present in dark sectors featuring color-charged particles and lead to a substantial impact on the theoretically predicted DM abundance. Therefore, similar consequences are expected for comparable scenarios in relevant theoretical or experimental analyses, which we would like to raise awareness for.
The strength of the potential acting on the initial/final state particles depends on the color configuration (i, j)/(i , j ) of the product representation according to Eq. (3.4). Thus, for an unambiguous initial and final state potential, the amplitude has to be decomposed into its irreducible representations. For the two relevant processes discussed in this article, we find M = P [1] M + P Using the fact that , we find for the squared matrix element of the process (X + X) [3] → B(XX) [3] + g [3] P [3] ijmn P [3] i j m n [M trans ] a i,i ,j,j Since we only consider bound state formed from particles equal in mass, we have η 1 = η 2 , which in turn implies that the squared amplitude for the BSF (X +X) [3] → B(XX) [3] +g [3] vanishes. Note that this result is independent from the approximation M 1 = M 2 , since any contribution involving M 1 vanishes individually. The remaining BSF process (X + X) [6] → B(XX) [3] + g [6] yields P Nevertheless, as discussed in [34], α NA s and α B g typically differ only up to O (10%), resulting in a suppression of this contribution to the BSF cross section by a factor of O (100). As a consequence, we do not include particle-particle bound states in our analysis.
The same conclusion can be derived by considering3 ⊗3 initial states.

B Boltzmann Equations including Bound State Decays via Constituent Decays
sector number density where {i, j} ∈ χ, X, X † and while all other channels are unaffected, i.e. σ ij v ij eff = σ ij v ij , for (i, j) = (XX † ), (X † X) .
We are now in the position to investigate the impact of a non-zero DM-SM coupling g DM , which allows for Γ B→Xχq = 0. On the one hand, Γ B→Xχq = 0 gives an additional positive contribution to the effective annihilation cross section of DM in Eq. (B.7). On the other hand, Γ B→Xχq = 0 increases the total bound state interaction rate Γ B,tot , thereby decreasing the effective annihilation cross section of DM. In order to analyze those two competing effects, we inspect the variation of the bound state contribution to the effective annihilation cross section with Γ B→Xχq Here, we take into account BSF+Sommerfeld processes and, in comparison with Fig. 5, also include effects of Γ B→Xχq = 0, more precisely we use Γ B→Xχq 2Γ X→χq . The dashed white lines are in correspondence of the region where 10 −7 g DM 10 −2 . Below this interval, it is not possible to ensure chemical equilibrium between the unbound particles in the dark sector and the co-annihilation freeze-out assumption breaks down. ∆. As a consequence, σ BSF v eff (g DM = 0) = 0 provides the smallest possible bound state contribution to the effective annihilation cross section of DM.
In this sense, the results presented in this paper, assuming a negligible three-body decay width of the bound state, provide the most conservative estimate of the effects of bound states in the context of the DM relic density. Furthermore, we verified that even a sizable three-body decay width does not alter our results presented in section 7. In order to do so, we approximate the three-body decay width by the sum of the decay width of the constituent particles Γ B→Xχq 2Γ X→χq . We present the results for the coupling g DM that leads to the observed relic density in Fig. 14.
In comparison with Fig. 5, which presents the same results assuming Γ B→Xχq = 0, we only find sizable effects for small mass splittings ∆m 30 GeV and DM masses above a few TeV. These regions are already clearly excluded by direct detection experiments as can be seen from Figs. 10, 11 and 12.