Cornering Colored Coannihilation

In thermal dark matter models, allowing the dark matter candidate to coannihilate with another particle can considerably loosen the relic density constraints on the dark matter mass. In particular, introducing a single strongly interacting coannihilation partner in a dark matter model can bring the upper bound on the dark sector energy scale from a few TeV up to about 10 TeV. While these energies are outside the LHC reach, a large part of the parameter space for such coannihilating models can be explored by future hadron colliders. In this context, it is essential to determine whether the current bounds on dark matter simplified models also hold in non-minimal scenarios. In this paper, we study extended models that include multiple coannihilation partners. We show that the relic density bounds on the dark matter mass in these scenarios are stronger than for the minimal models in most of the parameter space and that weakening these bounds requires sizable interactions between the different species of coannihilation partners. Furthermore, we discuss how these new interactions as well as the additional particles in the models can lead to stronger collider bounds, notably in jets plus missing transverse energy searches. This study serves as a vital ingredient towards the determination of the highest possible energy scale for thermal dark matter models.


Introduction
The thermal hypothesis, or the assumption that the dark matter used to be in thermal and chemical equilibrium with the Standard Model (SM) in the early Universe, tightly links the dark matter relic density to the strength of its interactions with the SM particles. The associated models hence predict a plethora of experimental signatures from colliders, direct, and indirect detectors. For minimal models, where the dark matter is the only new particle, the relic density bounds on the masses of the dark sector particles typically lie around a few TeV and these theories are therefore under siege at the current experiments [1][2][3][4]. One of the most efficient ways to loosen the existing limits is to introduce a new dark sector particle that will accelerate the dark matter depletion through the so-called coannihilation mechanism [5]. If this new particle is strongly interacting, coannihilation would allow the dark matter to reach masses of up to 10 TeV without overclosing the Universe [6][7][8][9][10][11][12][13][14][15]. This scenario is however significantly constrained by the LHC, which can probe O(10%) relative mass splittings between the dark matter and its coannihilation partner. While small mass splittings are particularly challenging to explore, a future 100 TeV collider should be able to probe most of the remaining regions of the parameter space [10,[16][17][18]. Although these projections seem extremely encouraging, it is essential to keep in mind that they have been derived from simplified models involving only the dark matter and a single coannihilation partner. It is therefore crucial to determine whether the current limits on the masses of the dark sector particles still apply to more complex scenarios.
A simplified model of dark matter coannihilation can be extended by adding either new mediators between the SM and the dark sector, new dark matter candidates, or new coannihilation partners. The first possibility is a straightforward way to loosen the constraints on the dark matter mass by adding new annihilation channels without increasing the new physics couplings. For these extended models, the constraints from perturbative unitarity, colliders, or relic density can be considerably weaker than the ones derived using simplified models, and are highly model-dependent [19,20]. We therefore do not study this configuration here. Conversely, as we will show in the rest of this paper, meaningful and generic constraints can still be derived for models with either multiple dark matter candidates or multiple coannihilation partners [21][22][23][24]. Although both scenarios can be studied using similar approaches, the parameter space for models with multiple coannihilation partners is larger since there are fewer restrictions on the partners' quantum numbers. In what follows, we will therefore focus on these types of models, keeping in mind that our techniques can be straightforwardly applied to models with multiple dark matter candidates.
In this study, we focus on models involving one dark matter candidate (DM) and an arbitrary number of strongly interacting coannihilation partners X i , close in mass to the dark matter. In this type of models, X i self-annihilates through strong interactions, thus causing an efficient depletion of the dark matter by shifting the DM-X i equilibrium. As shown in [6,9,11,12], the bounds on the dark matter mass from these models are among the weakest for thermal dark matter models, and can reach up to O(10) TeV. In the rest of this study, we build on the framework introduced in [6], where we assume that the dark matter relic density is entirely determined by the self-annihilation rate between the X i and all the processes involving the dark matter can be neglected. Note that in such scenarios the direct and indirect detection signals are suppressed and do not need to be considered. The dominant model-independent constraints on these type of models will therefore be the relic density requirement and the collider bounds on the X i .
Requiring the dark matter relic density to agree with the current observations [25] forces the dark matter mass into a narrow band whose central value highly depends on the mass splitting between the X i and the dark matter. For simplified models with a single coannihilation partner, the allowed dark matter mass ranges from O(100) GeV to about 10 TeV for low DM-X i splittings. Since the dark matter relic density can be increased ad libitum by introducing additional novel dark matter candidates, the lower bound on the dark matter mass set by relic density can be easily relaxed. The upper bound, however, is much more robust and can therefore be used to determine the range of energies that should be explored by the future experiments. Here, we show that the upper bounds on the dark matter mass derived for minimal coannihilating models are still valid when introducing new coannihilation partners, unless the "mixed" annihilation rates between different species of X i are particularly large. Moreover, we find that, in order to test the robustness of the simplified model constraints against introducing any number of X i , it is sufficient to consider models with only two coannihilation partners. This result considerably simplifies the study of non-minimal thermal dark matter scenarios.
In what follows, we derive the upper bounds on the mixed annihilation rates X i X j → SM SM up to which the simplified model constraints still apply. We study a representative set of models where the dark matter coannihilates with two strongly interacting particles X 1 and X 2 . We model the X 1 X 2 → SM SM process using only renormalizable interactions, that can be mediated by either a SM particle, a new particle in the s-channel, or a new particle in the t-channel. This approach allows us to express the upper bounds on this mixed annihilation rate in terms of constraints on the couplings between the different X i and the masses of the new mediators. In particular, we characterize several specific kinematic configurations where introducing new coannihilation partners can considerably enhance the dark matter depletion rate even for moderate couplings. We finally explore how introducing the new X i and their associated vertices affects the collider searches at the LHC and at a future 100 TeV accelerator.
The applications of our results are manifold. First and foremost, demonstrating the robustness of the current simplified model constraints is a necessary condition for an extensive search program for coannihilating dark matter at current and future hadron colliders. On the theoretical side, exploring models with extended dark sectors is significantly simplified as, in most of the parameter space, only one or two coannihilation partners need to be considered at the same time. In particular, our results are directly applicable to SUSY models where the gluino and/or the different squark flavors coannihilate with the neutralino. Note that, for squark-neutralino coannihilation, since flavor constraints strongly limit the annihilation rates between different sfermion species, the existence of multiple degenerate squark or slepton flavors can lead to particularly tight relic density bounds [26][27][28]. Another major class of models with multiple coannihilation partners is Kaluza-Klein theories where the dark matter is the lightest Kaluza-Klein particle and can coannihilate with higher modes that are nearly degenerate with each other [29].
Our study of models with multiple coannihilation partners is outlined as follows. In section 2, we discuss how the relic density changes when we increase the number of coannihilation partners X i in a given model, and we derive a set of upper limits on the mixed X i X j interactions. In the same section, we also introduce the simplified models of coannihilation that we are going to use in the rest of this work. We then describe all the possible tree-level mixed annihilation processes in more detail in section 3 and identify the regions of the parameter space for which adding new coannihilation partners can increase the dark matter depletion rate. Informed by these results, we discuss the collider constraints for these models in section 4. Finally we conclude in section 5.

Relic density for multiple coannihilation partners
In this section we discuss how increasing the number of coannihilation partners affects the dark matter relic density. In particular, we derive the conditions under which adding new coannihilation partners significantly increases the dark matter annihilation rate. Since our goal is to determine how heavy the dark matter mass can be in thermal scenarios, we only consider strongly interacting coannihilation partners, that typically lead to the largest annihilation rates in the dark sector. We follow the approach described in [6], focusing on several representative models that allow to derive generic conclusions about strongly interacting coannihilating dark sectors. In what follows, we first discuss the dependence of the effective dark matter annihilation cross-section on the annihilation rates of the additional partners and what ingredients are necessary for these partners to loosen the relic density bounds on the dark matter mass in any given model. We then describe the simplified models that we will use throughout this paper and present a few examples of how introducing additional coannihilation partners can affect the dark matter relic density constraints.

Multiple partners and dark matter annihilation
Here, we study how introducing multiple coannihilation partners modifies the dark matter effective annihilation rate in generic coannihilating models. We consider a scenario with one dark matter candidate with g DM degrees of freedom and N coannihilation partners X i with g i degrees of freedom. We denote the relative mass splittings between the X i and the dark matter by ∆ i = Neglecting the dark matter self-annihilation rate, the effective annihilation rate of the dark matter [5,30] in the non-relativistic approximation is where σ ij is the X i X j annihilation rate and g eff and K i are defined by In order to acquire an intuition on how the number of coannihilation partners affects the effective rate from equation (2.1), we first assume that all the X i are degenerate. We denote the mass splitting between the DM and all the X i as ∆. If only X i X i self-annihilation, with cross-section σ XX , is allowed, equation (2.1) becomes When ∆ approaches zero, which usually corresponds to the largest effective annihilation rates, this equation simplifies to which decreases as 1/N in the large N or small g DM limit. In the absence of mixed interactions between different species of dark sector particles, introducing additional copies of the coannihilation partner in any given model thus tightens the relic density constraints at low ∆. This rather counter-intuitive behavior has already been pointed out in the context of flavor violation in the squark sector [31].
At large ∆, the dynamics of coannihilation drastically changes. Although increasing the number of coannihilation partners still makes X i X i self-annihilation more difficult, the limiting process is now the conversion of DM into X i . Since adding new coannihilation partners increases the number of possible final states for this process, the effective dark matter annihilation rate will also increase. This phenomenon can be observed analytically by taking the large ∆ limit in equation (2.3), which gives In this regime, the increase of the effective annihilation rate is thus linear with the number of coannihilation partners.
Let us now introduce mixed X i X j → SM SM annihilation processes with cross-sections all equal to a given σ mix . Now, equation (2.3) becomes Since allowing for mixed annihilation does not modify the DM-X conversion rate, σ eff still increases with N at large ∆. In the small ∆ limit, however, the total annihilation rate is now nearly independent of N . In particular, if the X i annihilate with each other indifferently (σ mix = σ XX ), at ∆ = 0, the effective dark matter annihilation cross-section is the same as in a model with only one coannihilation partner. Conversely, when the mixed annihilation dominates over the self-annihilation (σ mix > σ XX ), we observe an increase of σ eff compared to a single-partner model. Interestingly, as long as there is more than one coannihilation partner in the model, this increase will depend only weakly on the actual value of N . This behavior indicates that focusing on models with only two coannihilation partners could be an efficient way to estimate how much σ eff can increase in generic coannihilating models.
As found in [6], the current LHC bounds on m DM and ∆ combined with the relic density constraints exclude values of ∆ down to about 10% for simplified models with a single coannihilation partner. Understanding the behavior of the dark matter annihilation rate in this low ∆ region is therefore crucial to inform the future collider search program. In the rest of this work, we thus focus on the ∆ ≈ 0 region and derive a set of sufficient conditions for the simplified model constraints on the dark matter mass to still hold in scenarios with multiple coannihilation partners X i . We will briefly comment on the effect of a non-zero ∆ on models with two different coannihilation partners in section 2.2. As already found in the simple case study above, we establish that the only way to invalidate the current simplified model constraints at low ∆ is to introduce sizable mixed annihilation rates between the X i . We will thus derive a set of upper bounds on these mixed rates for three different scenarios: N identical coannihilation partners, two different coannihilation partners, and finally N different X i . The latter is the completely general case and hence serves as our main result.

Multiple equal species
We first consider a simple scenario similar to the one discussed above, with one dark matter candidate that does not self-annihilate, and N identical coannihilation partners X i with the same numbers of degrees of freedom g X and self-annihilation cross-sections σ X . The mixed annihilation cross-sections between the different X i are assumed to be all identical and equal to σ mix XX . The effective dark matter annihilation rate in the small ∆ limit then reads This rate needs to be compared to the effective annihilation rate of a model with only one species, which in this case is defined by In order for σ eff N to be smaller than σ eff X for N ≥ 2 we therefore need In our models, the highest number of degrees of freedom for the dark matter is g DM = 6 for a complex vector boson, while the lowest possible number of degrees of freedom for a strongly interacting X is g X = 6 for a complex scalar triplet. The upper bound on the mixed annihilation rate will then range from σ X 8 for a complex vector dark matter candidate, a complex scalar triplet X, and N = 2, up to values close to σ X for large N and g X g DM . Since the X i can self-annihilate through the strong interaction, obtaining a large enough σ mix XX to break condition (2.9) requires either particularly large couplings or specific kinematic features such as resonances, interferences, or a suppressed (e.g. p-wave) self-annihilation cross-section.

Two different species
We now consider a model with one dark matter candidate and two coannihilation partners X 1 and X 2 that can have different properties. The effective annihilation cross-section for this model is This cross-section needs to be smaller than the single effective annihilation cross-section for a single-partner model with either DM and X 1 or DM and X 2 , given by Now we can freely assume σ eff X 1 > σ eff X 2 , which leads to a condition of the form Note that this equation reduces to the constraint in equation (2.9) for N = 2 and equal X 1 and X 2 , as expected.

Multiple different species
We finally consider the most general scenario, with one dark matter candidate and N coannihilating partners that are allowed to have different properties from each other. We can rewrite equation (2.1) as (2.13) Here, the second term sums over all mixed annihilation cross-sections between the coannihilation partners X i and X j . This time, σ eff X 1 ···X N may not exceed max σ eff X 1 , . . . , σ eff X N , with σ eff X i the effective annihilation cross-section for a model with only the dark matter and X i defined in (2.11). Although this constraint naively translates into heavily modeldependent bounds on the different mixed rates, it is also automatically satisfied when all the submodels of the form DM + X i + X j satisfy equation (2.12). The detailed derivation of this property is shown in Appendix A. As a consequence, all the possible mechanisms that could loosen the existing simplified model bounds in models with multiple coannihilation partners already appear in models with two X particles. This is our main result and considerably simplifies our analysis, since it allows us to focus exclusively on such models for the rest of this work.
We now apply the requirements derived here to specific dark matter models. The formalism used throughout this paper, in particular the different models that we are going to study, is detailed in the next section. We highlight the importance of the mixed X i X j interactions through several examples and briefly comment on the cases where the mass splittings ∆ i are different from each other.

Models of the coannihilating dark sector
In this work we extend the scope of the previous colored dark sector studies [6,9,[11][12][13] to include multiple coannihilation partners. We follow the methodology described in [6] and consider simplified models where the dark matter is a pure SM singlet that does not self-annihilate. The dark matter is in thermal and chemical equilibrium with at least two coannihilation partners X i that are strongly interacting. As shown in [6], in such models, the self-annihilation of the X i through the strong interaction significantly contribute to the dark matter depletion. These processes alone can in fact loosen the upper bound on the dark matter mass from a few TeV up to more than 10 TeV. In this study, in addition to the SM couplings, we will also introduce interactions that allow for mixed annihilation channels of the form X i X j → SM SM. We however neglect the influence of the X i DM → SM SM coannihilation channels (see reference [7]) on the final relic density. Since these processes are still necessary to ensure that the X i decay and are in equilibrium with the dark matter, we introduce them as effective operators of the form DM X i SM SM. The final Lagrangian for a given model will then be of the form (2.14) The kinetic and mass terms of X and DM, L X and L DM as well as the effective operator describing the DM-X i interactions are taken from [6]. In what follows, we will allow the X i to be either complex scalars or Dirac fermions, in the triplet, sextet or octet representation of SU (3). The corresponding L X i for X being either a scalar S or a fermion ψ are of the form where i, j are color indices and the T a R matrices are the generators for the color representation R of X. The covariant derivatives are given by We assume that the dark matter field as well as X i are protected by a Z 2 symmetry. Since the dark matter quantum numbers will only enter in our final results through the effective number of degrees of freedom, we put no restriction on the spin of the DM.
While, as in [6], the DM-X i interactions are modeled by a suppressed effective operator, this formalism is not suitable for the X i X j SM SM interactions since in our study the associated annihilation rates can be large. In what follows, we will therefore fully resolve these mixed interactions by introducing a mediator that can be either a SM particle, a new physics s-channel, or a new physics t-channel particle. In the last two cases, we allow the mediator to be either a scalar, a fermion, or a vector, in the triplet, sextet, or octet representation of SU (3), and choose the color and the spin that leads to the largest dark matter effective annihilation cross-section for each model. We introduce the L X i +X j associated to the three possible mediator configurations in section 3.1. We provide all the models discussed in this work and in [6], including the mixed interactions, in the updated version of our FeynRules v2.3.24 [32,33] model package [34]. Using a Mathematica notebook each of the specific models can be generated in both UFO [35] and CalcHEP v3.6.25 [36] format and then be further used to do collider studies or calculate the thermal relic abundance.
In figure 1, we show the effect of additional coannihilation partners, with and without mixed interactions, on the relic density constraints for a few of the models described in equation (2.15). Here, we start from a model involving one fermionic dark matter candidate and one fermionic color octet X F8 (blue band) and add new coannihilation partners either neglecting mixed interactions (left figure) or saturating condition (2.12) (right figure). For all models we compute the relic density using micrOMEGAs v4.3.2 [37], directly inputting the cross-sections for the mixed processes into the code without assuming a specific type of interaction. As discussed in section 2.1, introducing extra particles, whether with the same or different quantum numbers always lowers the allowed dark matter masses when ∆ is small and the mixed condition (2.12) is satisfied. However, for ∆ 7% in the absence of mixed interactions there is a crossover and introducing new coannihilation partners increases the allowed dark matter mass, as expected from equation (2.5). Finally, when saturating the condition (2.12), the allowed dark matter masses become independent on the number of coannihilation partners at ∆ ≈ 0 and increase with the number of X particles for larger ∆. Importantly, for any model the largest allowed dark matter masses are still obtained at ∆ = 0.
For the former analysis as well as for the rest of this paper, we restrict our study to the case of ∆ 1 ∆ 2 . In order to illustrate how relaxing this assumption changes the bounds on the DM masses, we compute the relic density constraints on a model with two fermionic color octet coannihilation partners, X F8 and X F8 , with mass splittings ∆ 1 and ∆ 2 with the dark matter. The dark matter masses leading to the correct relic density are shown in figure 2  Figure 1. Relic density constraints in the ∆ versus m DM plane for models where DM is a Dirac fermion. We show the parameter space regions that agree with the relic abundance measured by Planck. We compare models with multiple coannihilation partners against the benchmark model with a single partner X F8 (solid blue line). We either increase the multiplicity of X F8 or introduce an additional partner X C3 . The left figure shows the contours in the absence of mixed annihilation channels, while the right figure shows the models with the mixed annihilation cross-section saturating condition (2.12).
as a function of ∆ 1 and ∆ 2 . In this figure the straight horizontal and vertical lines indicate the regions where one of the coannihilation partners decouples. In the absence of any X 1 X 2 annihilation process and for ∆ 1,2 10%, bringing the ∆ i closer to each other leads to a decrease in the DM mass, and thus to tighter relic density constraints. For ∆ 1,2 10%, the opposite behavior occurs. These results are a direct generalization of the behavior observed on the left panel of figure 1 since the ∆ i ∆ j and the ∆ 1 = ∆ 2 cases correspond to the DM F + 1X F8 and DM F + 2X F8 models respectively. A similar reasoning can be applied to the scenario where the X 1 X 2 annihilation rate saturates the condition (2.12), shown on the right panel of figure 2. Here, by construction, for low values of ∆ 1,2 the allowed dark matter mass for ∆ 1 = ∆ 2 is similar to the one obtained when one of the coannihilation partners decouples. For larger ∆ i , however, an increase in the allowed dark matter mass of up to about 100 GeV can be observed in the ∆ 1 ≈ ∆ 2 region. The diagonal symmetry observed for both the left and right panels of figure 2 stems from the fact that the two coannihilation partners in our model have identical quantum numbers. For models where this is not the case, the diagonal would be shifted but the general features of figure 2 are conserved. The results obtained throughout this section for identical ∆ i are therefore directly applicable to scenarios with different ∆ i and we can safely assume identical mass splittings for the rest of this study.
Note that, since we present figures 1 and 2 for illustrative purposes, we computed only perturbative annihilation rates. In principle, these rates get modified by non-perturbative correction factors from the Sommerfeld effect [38,39] and bound state formation processes [12,40]. As discussed in detail for models with a single coannihilation partner in [6], these effects generically lead to an increase in the dark matter mass allowed by the relic density requirement. This shift in the dark matter mass, however, will occur for both the single and multiple partner models so the relative positions of the relic density bands shown in figure 1 will not be significantly modified. Moreover, the results derived in section 2.1 do not depend on how the annihilation cross-sections are computed and are thus valid whether or not non-perturbative effects are taken into account. We confirm this result through a quantitative study in section 3.2.
As discussed in section 2.1.3, in order to study the validity of the usual simplified model constraints for models with multiple coannihilation partners, it is sufficient to consider models with only two X particles. In the next section we will consider all possible models with one dark matter candidate and two coannihilation partners X 1 and X 2 that are charged under SU (3). Besides the gauge interactions, we introduce additional vertices and particle content in order to allow for the mixed X 1 X 2 → SM SM interaction. The corresponding models will therefore be characterized by the quantum numbers of X 1 and X 2 , the dark matter mass m DM , the relative mass splittings ∆ 1 and ∆ 2 between the dark matter and X 1 , X 2 , and the parameters associated with the mixed interactions. In what follows, as discussed in section 2.1, since we are interested in the region of parameter space giving the loosest limits on the dark matter mass, we assume that ∆ 1 = ∆ 2 = 0. As derived in section 2.1, the maximum value on the mixed annihilation rate up to which the simplified model constraints are valid is then This condition is the main result of our work and serves to ensure that at low ∆ i the introduction of additional coannihilation partners does not increase the maximal dark matter mass allowed by the relic density constraint.
In what follows, we fully resolve the structure of the X 1 X 2 → SM SM interaction and investigate the impact of condition (2.17) on the corresponding parameters. Since, in order to saturate this condition, the mixed annihilation rate needs to be large, we will consider only tree-level and renormalizable interactions. This restriction leaves us with three possible topologies: annihilation via a SM mediator, s-channel annihilation via a new physics mediator, and t-channel annihilation via a new physics mediator. The mixed interaction will therefore be characterized by at most two coupling constants as well as the mass and the quantum numbers of the mediator. A detailed study of these three scenarios and the associated constraints is presented in the next section.

Mixed annihilation
In the previous section, we derived a sufficient condition for the constraints from simplified models with a single coannihilation partner to still apply to models with multiple coannihilation partners. This condition can be written as an upper bound on the annihilation rates for "mixed" processes of the form X i X j → SM SM. Moreover, we showed that to understand whether the coannihilating simplified model constraints still apply to a model with an arbitrary number of coannihilation partners, it is sufficient to apply this requirement to each X i X j pair individually. In what follows, we explore how this condition on the mixed annihilation rates constrains the masses and couplings in the dark sector for models with two coannihilation partners. To this end, we first construct all the possible tree-level renormalizable interactions that can lead to mixed coannihilation processes, introducing new particles and new couplings if necessary. For each of these processes and for representative choices of quantum numbers for the X i and the dark matter, we then derive the bounds on these new masses and couplings within which the coannihilating simplified model limits on the dark matter mass are still valid.

Modelling the mixed interaction
Adding coannihilation partners to a given model opens new annihilation channels that will modify the dark matter depletion rate. In particular, aside from the self-annihilation channels X i X i → SM SM studied in [6], processes of the form X i X j → SM SM can now take place. As outlined in section 2.1, in order to establish whether introducing additional coannihilation partners can significantly loosen the bounds on the dark matter mass, it is essential to model these mixed processes and estimate their strength.
In the rest of this section, we focus on all the possible X i X j interactions that could significantly modify the dark matter annihilation cross-section. As discussed in section 2.1, we consider models with two coannihilation partners X 1 and X 2 and choose both of these partners to be strongly interacting since we focus on the models typically associated to the loosest relic density bounds. Following the approach described in [6], we assume that the thermal and chemical equilibrium between a given X i and the dark matter is ensured by an effective operator of the form X i DM SM SM, with the DM being an SM singlet. In the rest of this work, we will assume that this operator is suppressed and can be neglected compared to the other annihilation channels. The Lagrangian for a given model can then be written as where L DM , L X 1 , L X 2 , L DM−X are given in section 2.2. We have implemented several examples for L X 1 −X 2 in the Mathematica package [34] introduced in [6] for single partner models. As seen in section 2.1, the rates of the mixed and self-annihilation processes need to be comparable to lead to an increase of the DM annihilation cross-section. In what follows, we therefore ignore possible non-renormalizable or loop-level interactions between X 1 and X 2 and focus on the dimension four tree-level interactions susceptible to lead to X 1 X 2 → SM SM processes. These interactions are of three types: X 1 X 2 SM, X 1 X 2 M s with M s being a new physics s-channel mediator, and X 1,2 M t SM with M t being a t-channel mediator interacting with both X 1 and X 2 . In the rest of this work, we study these three possible scenarios separately for various color representations of X 1 and X 2 . Since our goal is to establish what the loosest constraints on the dark sector can be in these models, we will consider only the Standard Model particles giving the weakest LHC bounds, i.e. quarks and gluons. We would like to emphasize, however, that the procedure outlined in this section is generic and is expected to lead to similar results for a wide range of models involving weakly interacting particles.
For two-to-two annihilation processes into quark or gluons, the allowed color representations for the new physics particles naively range from 1 to 27. In fact, since two different particles cannot annihilate into two gluons at tree-level, the only allowed SU (3) representations for X 1 and X 2 are 1, 3, 6 and 8. In what follows, we generally choose the X partners to be either triplets or octets of SU (3) and the mediators to be either singlets, triplets, or octets. Although sextets are commonly encountered in new physics models [41,42], for most processes, their associated color factors are similar to the ones for octets. We therefore consider color sextets only for processes where color octets are forbidden by gauge invariance, such as s-channel annihilations to q q orqq. We do not impose any particular restriction on the spin of the two X particles, that can be either scalars, fermions, or vectors. For the sake of clarity, we do not consider all the possible spin configurations for the new physics particles and, for a given mixed process, just choose the spin and coupling structure that give the least suppressed amplitude. Finally, note that, in our framework, the dark matter does not self-annihilate and its quantum numbers contribute to the effective annihilation cross-section only through the number of degrees of freedom g DM . Since according to equation (2.17) a larger g DM enhances the contribution of the additional X j to the effective annihilation cross-section, we choose DM to be a complex vector, which corresponds to g DM = 6.
Although we neglect processes involving weak interactions, it is important to notice that specific UV completions of the X 1 -X 2 -SM-SM interaction can lead to constraints on the SU (2) quantum numbers of X 1 and X 2 . These constraints notably arise when X 1 and X 2 annihilate to qq via an s-channel mediator M s since the SU (2) quantum numbers of the quarks depend on their chiralities. After EWSB, scenarios where the coannihilation partners are SU (2) multiplets can involve a large number of coannihilation partners and channels and need to be treated with care [2,43,44]. Since the dark sector particles in our models typically have multi-TeV masses, we can however safely treat all the components of a given SU (2) multiplet as identical copies of the same particle that -since we neglect the effects of the weak interaction -do not mix with each other. In this approximation, according to the results from 2.1.3, scenarios where X 1 and X 2 are SU (2) multiplets can be studied by considering only a single component of each multiplet. Hence, in what follows, in models where X 1 and X 2 are required to be charged under SU (2), we will only study the annihilation of the electrically neutral components of X 1 and X 2 through the neutral component of the mediator, effectively treating these particles as SU (2) singlets.
We now introduce the different types of X 1 X 2 → SM SM processes, classifying them by their mediator: SM particle, new physics particle in the s-channel, or new physics particle in the t-channel. For each category, we present the Lagrangians for the models we study and discuss in detail all the associated annihilation processes.

SM mediator
A common new physics scenario is the existence of a coupling between two dark sector particles, X 1 and X 2 , and a Standard Model field. For the models studied here, since we are considering only annihilations into quarks or gluons, this Standard Model field can only be a quark. The SU (3) representations for X 1 and X 2 allowed by gauge invariance are therefore (3,3), (3,6), (3,8) and (6,8) as well as their conjugates. Similarly, the only allowed spin configurations are the ones where one of the X is bosonic and the other one is fermionic. For these choices of quantum numbers, introducing the new X 1 -X 2 -quark interaction leads to the additional coannihilation processes shown in figure 3. Note that a new t-channel X i X i → SM SM annihilation diagrams appears, which will either increase or decrease the self-annihilation rate of the coannihilation partners, depending on the quantum numbers of the X i .
For this study, we choose X 1 to be a complex scalar and X 2 to be either a Dirac or a Majorana fermion. Note that, if X 2 is a Majorana fermion, either of the new annihilation channels X 2 X 2 → q q,qq opens. The existence of this new channel leads to an increase in the X 2 self-annihilation rate that makes is more difficult for the mixed annihilation rate to dominate compared to the Dirac fermion scenario. We focus on a few representative Figure 3. The first three diagrams show the mixed annihilation processes when a coupling between the two X's and a SM quark is present. This coupling also introduces a new contribution to the self-annihilation rates of X 1 and X 2 as shown in the fourth diagram. models for which the X 1 -X 2 interaction Lagrangian is shown below In these Lagrangians the K u ij and their properties are defined in [45], and the color factor K aui for the 3 − 6 − 8 interaction is where j, k, and l run from one to three and i, u, a are the indices of the external triplet, sextet, and octet particles respectively. Here, q can be either an up or a down-type quark depending on the electromagnetic charges of X 1 and X 2 . Besides the particle masses, the only new physics parameter for these models is the coupling y NP , which will typically be of order one as all particles involved are colored.

New physics mediator (s-channel)
Instead of directly coupling to a SM particle, X 1 and X 2 can also interact with the SM through a new physics s-channel mediator M s . The corresponding interaction requires introducing two new trilinear couplings, for the M s -X 1 -X 2 and the M s -SM-SM vertices, thus constraining M s to be Z 2 even. An interesting aspect of this model is that the new interactions open the mixed annihilation channel shown in figure 4 without introducing any new self-annihilation diagrams. Hence, such s-channel models could potentially lead to significant increases in the total dark matter annihilation cross-section with respect to models with a single partner, especially near the resonant m M ≈ m X 1 + m X 2 region. The increase of the annihilation rate near the resonance, however, is primarily due to the introduction of M s in the model and could therefore also occur in a scenario with only one coannihilation partner that self-annihilates through a new physics mediator. In order to disentangle the effects of adding a new coannihilation partner from the effect of having a new mediator, our two-partner model therefore needs to be compared to a one-partner model with a new mediator that has similar properties as M s .
Since, for a new physics mediator, M s -SM-SM interactions involving gluons are forbidden at tree-level, we focus only on scenarios where X 1 and X 2 annihilate into quarks. These annihilation processes can occur if the mediator is either a scalar or a vector, and a singlet, triplet, sextet or octet of SU (3). Here, as mentioned at the beginning of this section, we consider only triplet and octet mediators and study models for which the mixed annihilation cross-section has a sizable s-wave component. For each X 1 X 2 pair, we choose the spin of the mediator that leads to the largest X 1 -X 2 annihilation cross-section. We therefore always use vector mediators for models with fermionic X 1 and X 2 and scalar mediators for models where the coannihilation partners are scalar fields. Mediators in the self-conjugate 1 and 8 representations of SU (3) are taken to be real while triplet mediators are constrained to be complex. These choices lead us to select the following representative models Additionally, for scalar X 1 and X 2 , all the M s -X 1 -X 2 interactions are parameterized by a dimensionful trilinear coupling that we write as A X = y X m S . This notation highlights the fact that perturbative unitarity bounds prevent the trilinear coupling to be larger than a few times the mass of the mediator [46,47].
As mentioned at the beginning of this section, in order to isolate the effects of introducing a new coannihilation partner for a given s-channel model, we need to compare this model to a scenario with only one coannihilation partner X that annihilates to qq via a mediator M s . In order for the two models to be comparable, M s needs to be as similar as possible to the mediator of the mixed coannihilation process M s . For models where this process is mediated by either a singlet or an octet, the quantum numbers of M s and M s can be taken to be exactly identical. For models where the X 1 X 2 annihilation is mediated by a triplet, however, M s cannot be reused to mediate the XX self-annihilation due to gauge invariance. In these cases, we consider both the cases of a singlet and an octet M s , that has the same spin as M s . In order to further ensure a fair comparison between the twoand one-partner scenarios, the couplings y X and y SM of M s to XX and the SM respectively as well as the spin structure of the corresponding vertices are taken to be the same as for their mixed counterparts. For example, a model where X 1 and X 2 are both scalar triplets can be compared to a model with only one scalar triplet X and a mediator M C1 , where Figure 4. The mixed annihilation diagram when a new physics mediator is present. The left diagram describes the case of an s-channel mediator whereas the middle diagram describes t-channel mediators. The right diagram denotes the modification of the self-annihilation rates of the coannihilation partners for a t-channel mediator.
the new physics interactions are parameterized by the Lagrangian L X C X C M C1 shown in equation (3.4).

New physics mediator (t-channel)
In this scenario, two new vertices M t -X 1 -SM and M t -X 2 -SM need to be introduced in order for X 1 and X 2 to annihilate to SM particles in the t-channel. Here, M is required to be Z 2 odd and part of the dark sector. The mediator should therefore be heavier than both the dark matter and X for our choice of parameters. We will in particular focus on regions of parameter space where the splitting between M t and X 1 , X 2 is sufficiently large to avoid chemical or thermal equilibrium between M t and the dark matter before and during freeze-out.
The introduction of the mediator and X 1 , X 2 with their associated vertices opens the new annihilation channels shown in figure 4. Interestingly, in this class of models, it is impossible to increase the X 1 X 2 annihilation rate without also increasing the self-annihilation rates of X 1 and X 2 . Consequently, even for large couplings, any increase in the dark matter annihilation cross-section in the two-partner model compared to models with only either X 1 or X 2 is primarily due to the introduction of a new mediator. Since, as mentioned in the introduction, studying the effect of extra mediators on the dark matter depletion rate is beyond the scope of this work, we carefully factor out any mediator contribution to the dark matter annihilation rate in our subsequent study, as for the s-channel models studied in section 3.1.2. We describe the associated procedure in more details in section 3.4.
As in the previous scenarios, annihilation of X 1 and X 2 into two gluons is forbidden by gauge invariance and we focus on X 1 X 2 annihilation into either q q or q q. We therefore consider scenarios where either X 1 and X 2 are both complex scalars with a fermionic mediator, or X 1 and X 2 are both Dirac fermions with a complex scalar mediator. In this scenario, whether the fields are self-conjugate does not change the coupling strength for the interactions involved in the annihilation processes and we therefore made arbitrary choices. For all our models, we constrain X 1 , X 2 , and the mediator to be SU (2) singlets, their hypercharges being defined by gauge invariance for each model. Finally, as for the models with an s-channel mediator, we choose the new physics particles to be either color triplet or octets. The Lagrangians describing the X 1 and X 2 interactions with the SM for our choice of models are for scalar coannihilation partners and L X F3 +M S1 = y NP M S1 X F3,i q R,i + h.c.
for fermionic coannihilation partners. A model for the mixed annihilation is constructed by combining either two different Lagrangians with a same mediator or two copies of a same Lagrangian with different coupling strengths if X 1 and X 2 have the same quantum numbers.

Our procedure
In what follows, we describe how the effective annihilation cross-section changes when introducing new coannihilation partners for each of the models described above. More specifically, in addition to determining when condition (2.17) is verified, we also estimate how large the annihilation rate can become for models outside its range of validity. To this end, for each model with X 1 , X 2 , and possibly a mediator, we evaluate the ratio of the effective annihilation cross-section in the complete model over the value that this cross-section would take if either X 1 or X 2 is removed. Thus, in models where the X 1 -X 2 annihilation process requires a new physics mediator, any possible effect of this mediator on the self-annihilation cross-section of the X i will appear both in the numerator and the denominator. Our procedure therefore allows to factor out the effects of the mediator on the self-annihilation cross-section in the t-channel models studied in section 3.1.3.
For fixed particle momenta, the ratio of the effective annihilation rates, denoted by r, can be written as follows For a given choice of X 1 and X 2 , r depends on the masses and couplings of the new particles as well as on the quantum numbers of the dark matter and the mediator. Condition (2.17) now corresponds to r ≤ 1 and will translate into constraints on all these parameters. For models with more than two coannihilation partners, the range of validity of condition (2.17) is simply the intersection of the constraints corresponding to all the possible X i -X j combinations. The actual values of r, however, cannot be readily extrapolated from models with two coannihilation partners up to models with an arbitrary number of X i . Estimating these values for two-partner models nevertheless provides a resonable indication of how large annihilation rates can become in regions where (2.17) is violated.
Since r heavily depends on the momenta of the particles involved in the annihilation processes, in the rest of this study, we will consider a slightly modified version of this ratio in order to select the particle velocities most relevant around freeze-out times. More precisely, we average all the annihilation cross-sections over a Maxwell-Boltzmann velocity distribution, and define a new ratio R as follows with v being the velocity of an incoming particle in the center of mass frame and x ≡ m DM /T . Since our main goal is to investigate the relic density constraints on the dark sector, we evaluate R at x = m DM /T = 25, which is a typical value for the freeze-out temperature. As mentioned in section 2.1.1, since R increases with g DM , we choose the dark matter to be a complex vector in order to obtain the weakest possible relic density bounds. Now, naively, R is expected to only depend on the masses of the X i and the mediator, as well as on the new physics couplings. For ∆ 1 = ∆ 2 = 0, however, at perturbative level, R will depend on the masses of the new particles only through the ratio of the mediator mass m S over the mass m X of the X i . In the rest of this work, we will therefore study the dependence of R in m S /m X and in the new physics couplings. In scenarios where we include Sommerfeld corrections, since the value of α s in the QCD potential depends on the mass of X, we fix m X = 1 TeV and check that varying m X in the range 500 GeV ≤ m X ≤ 10 TeV does not quantitatively change our results.
In what follows, we will compute R for each of the models described in sections 3.1.1, 3.1.2, and 3.1.3 as a function of the new physics couplings as well as the ratio of the mediator mass over m X . We consider these three scenarios separately and, for each of them, discuss in detail the regions of parameter space where condition (2.17) is violated. For most scenarios, we consider only color triplets and octets and ignore non-perturbative effects. In order to assess the validity of these restrictions, we perform a complete study of the scenarios where the X 1 -X 2 annihilation is mediated by a SM quark, including both Sommerfeld corrections and models involving color sextets. We show that the Sommerfeld corrections do not qualitatively modify the constraints on the new physics parameters from condition 2.17. Moreover, for any given model involving a color sextet, for couplings of order α s and above, the values of R are always smaller than the ones obtained when the sextet is replaced by either an octet or a triplet, as long as these other color configurations are allowed. For the models involving a new physics mediator, these color configurations are usually permitted and these models will therefore not qualitatively change our results. We thus consider scenarios with color sextets only for the models with SM mediators shown in the next section.

Standard Model mediators
We first consider the case of mixed annihilation processes mediated by a SM quark. These processes require introducing a X 1 -X 2 -q interaction, that will give rise to three new mixed annihilation channels, depicted in figure 3. Interestingly, this new vertex also allows for new self-annihilation diagrams for X 1 and X 2 , as shown in figure 3. Since these new selfannihilation modes only appear when the mixed process is introduced, their impact on the self-annihilation cross-sections is not taken into account in condition (2.17). For the models studied in this section, we therefore slighly modify equation (2.17) to account for the new self-annihilation rates, as follows and the ratio R in equation (3.8) is modified accordingly. Here σ mod X i is the modified selfannihilation rate obtained when including the fourth diagram in figure 3.
α NP 3.8 α s α NP 4.9 α s X C3 X F8 6 32 α NP 5.3 α s α NP 7.0 α s X C6 X F3 12 12 α NP 2.1 α s α NP 3.9 α s X C6 X F8 12 32 α NP 3.7 α s α NP 3.9 α s X C8 X F3 16 12 α NP 2.6 α s α NP 6.3 α s X C8 X F6 16 24 α NP 6.9 α s α NP 7.6 α s Table 1. Values of the new physics coupling α NP that saturate the mixed condition (3.10) for new physics models where the mixed annihilation is mediated by an interaction with SM quarks. We marginalize over the other relevant parameters g DM , m X and v as described in the main text.
Since the models studied here do not include a new physics mediator, the ratio of the rates, R, only depends on the X 1 -X 2 -q coupling, α NP ≡ y 2 NP 4π . For each of the models described by Lagrangian (3.2) we compute R as a function of α NP and find the maximal value of α NP for which condition (3.10) is verified, or, equivalently, R ≤ 1. The results for the different models are summarized in table 1, both with and without Sommerfeld corrections. Details about our calculation of the Sommerfeld effect can be found in a previous work [39] and accompanying Mathematica package [48], as well as in appendix B.   For all the models listed in table 1, in order to significantly increase the dark matter depletion rate compared to single-partner models, we need either new physics couplings α NP much larger than the strong coupling or couplings of X 1 and X 2 to multiple quark flavors at a time. Since the second option is typically strongly constrained by flavor measurements, we conclude that for models where the mixed X 1 -X 2 annihilation occurs via a SM quark, increasing the number of coannihilation partners typically leads to stronger relic density bounds on the DM mass. In order to estimate the effect of α NP on the DM annihilation for theories with very large couplings, we also show R as a function of α NP in figure 5. We can see that, for large values of α NP , models involving color triplets are associated with the largest values of R, followed by models with sextets and octets. These results confirm our previous hypothesis about the dependence of R on the colors of the annihilating particles. Scenarios where mixed annihilations occur via new physics mediators also exhibit this feature, which is due to the fact that the X particles in the lowest SU (3) representations are also associated with the lowest annihilation rates and hence the tightest relic density bounds. Even for color triplets, however, adding new coannihilation partners loosens these bounds only for large couplings α NP α s . For such couplings, the pair-production of X at colliders can be significantly enhanced, leading to tighter collider bounds. We discuss this last point in section 4.

s-channel mediators
Here, we evaluate the contribution of the mixed annihilation rate to the effective DM annihilation cross-section for models where the X 1 X 2 → SM SM interaction occurs through a new physics s-channel mediator M s . We focus on a set of representative models, whose Lagrangians are shown in (3.4). As outlined in section 3.1.2, in order to factor out the effect of introducing the s-channel mediator M s , we compare the effective annihilation rates for a model with X 1 , X 2 , and M s to a model with either X 1 or X 2 and a mediator M s that couples to qq and X i X i . The mass of this mediator as well as the values of its couplings to the SM and the dark sector are taken to be the same as for M s . Besides the masses and degrees of freedom of the dark matter and the X partners, the DM depletion rate for these models is thus characterized by four parameters: m M , m X , α SM , and α X . Since, as shown in section 3.2, the Sommerfeld corrections should not significantly influence the final results, we consider only the perturbative (co)annihilation cross-sections, that depend on α SM , α X , and ρ = m M m X . In order to facilitate the interpretation of our results, we consider the ratioR(r, α max ), defined bỹ with R given in equation (3.8).
In order to understand how the mixed X 1 X 2 s-channel interactions modify the DM annihilation rate, we rewrite the annihilation cross-sections corresponding to the X i X j → M ( ) s → qq processes as where C ij is the color factor corresponding to the interaction, Γ ij is the total decay width of the mediator, S ij is a symmetry factor that is equal to 1 2 if exactly one of the initial states is self-conjugate and to one otherwise, andσ is the part of the cross-section that is the same for the one-and two-partner models. The decay width Γ ij depends on the new physics couplings in the following way where the Γ ij,SM and the Γ ij,X depend only on the spins and the color charges of the mediator and the X i , as well as on the masses of the dark sector particles. In what follows, since we are considering only perturbative theories, we assume that this width is narrow, enforcing the loose requirement that Γ ii , Γ ij < 1 4 m M , with Γ ii and Γ ij the decay widths of the mediator in the one-partner model leading to the largest annihilation rate and in the two-partner model respectively. Here, contrary to the SM mediator case studied in section 3.2, the new physics couplings appear both in the numerator and the denominator of σ ij . Consequently, just increasing these couplings does not necessarily lead to an increase in the cross-section for the mixed annihilation processes. This cross-section is in fact large only in the following regimes: • The resonance region m M ∼ 2m X : for non-relativistic dark sector particles, the annihilation cross-section reduces to For a reasonably narrow width Γ ij < 1 • The large coupling limit α X 1 with m M < 2m X : in this scenario, Γ X = 0 so the denominator of the cross-section does not depend on α X . σ ij can therefore be rewritten as and grows linearly with α X .
• The large coupling limit α X 1 (α SM 1) when Γ X (Γ SM ) is small: in this case, it is possible to show that, away from the resonance and in the narrow width approximation, the cross-section is maximized when the total width is equal to its maximal allowed value Γ = 1 4 m M . Optimizing over α SM , this value can be obtained for 14) and the cross-section is then
An equivalent formula can be obtained when optimizing over α X , by replacing α SM ↔ α X and Γ SM ↔ Γ X . The maximal coupling given in equation (3.14) can reach large values if Γ SM is small, which can in turn lead to large values for the cross-section. In this case, the latter would be only bounded by perturbativity.
In all these scenarios, for α X α s , the cross-sections for the new physics processes will dominate over the ones for the QCD processes. The ratio of the ratesR can then be approximated as This ratio reaches a maximum at the resonance, where it can be approximated bỹ (3. 16) Our estimates ofR in the regions where the new physics annihilation processes dominate thus show thatR is always less than one for the models where X 1 and X 2 have the same quantum numbers, as for these models, Γ ii = Γ 12 and C ii = C 12 andR can therefore be approximated asR (3.17)  Figure 6.R = 1 contours in the (α max , mM mX ) plane, withR defined in equation (3.11). The thick solid lines show the regions for which the width of the mediators are smaller than 25% of their mass. The particle content of each model is shown as a label next to the corresponding line and is of the form X 1 + X 2 + M 2-partner (+ M 1-partner ) where M 2-partner and M 1-partner are the mediators of the two-partner and one-partner models respectively. For the X F3 + X F3 + M W8 (+ M W8 ) model, the two-partner mediator mediates both the mixed and self annihilation of the X i . For models where X 1 and X 2 have different quantum numbers,R can be larger than one and is expected to be maximal in the resonant m M ≈ 2 m X region. We verify this hypothesis for the X F3 + X F8 + M W3 and the X C3 + X C8 + M C3 models, for which theR = 1 contours in the α max versus m M plane are shown in figure 6. Since we are exploring regions of the parameter space where the new physics couplings are large, the decay width of the mediators in these models can be comparable to their mass, which would indicate an underlying non-perturbative dynamics. In figure 6 we therefore also indicate the regions for which the decay width of all the mediators in our models are less than 25% of their mass. Applying this narrow width requirement constrains theR > 1 regions for all our models to be within about 30% of the resonance. Even when this criterium is relaxed, thẽ R > 1 domain remains relatively narrow, with the mediator mass remaining within 50% of the resonant region for α max < 2 α s .
In order to estimate how much our two-partner models can enhance the dark matter annihilation rate, we now fix α max so that it saturates the narrow width requirement Γ M < 1 4 m M , and plotR as a function of the mediator mass in figure 7. Here again,R sharply peaks near the resonance and its maximal values are close to ones predicted from equation (3.16). For most of our models, however,R remains less than two, which corresponds to an order one increase in the relic density bound on the dark matter mass. Significantly loosening Finally, we consider the possibility that, for models where X 1 and X 2 have the same quantum numbers, M s can couple not only to X 1 -X 2 but also to X 1 -X 1 and X 2 -X 2 . In this case, the total DM annihilation rate for the two-partner models can be further enhanced with respect to the one-partner model, even far from the resonant region. More specifically,R can be written asR where i can be either 1 or 2. For s-channel models, the QCD and the new physics processes interfere positively for the X i X i → qq annihilation so σ X i X i ≥ σ X 1 X 2 . Hence, we always haveR which always lead to values smaller than two for colored coannihilation partners. We confirm these results by showing theR = 1 contours for this model in figure 6 as well as the maximal values forR as a function of the mediator mass in figure 7. As shown in figure 6,R can now be larger than one further away from the resonance than in the X F3 + X F8 + M W3 and the X C3 + X C8 + M C3 models. Figure 7 shows however that the increase in the dark matter annihilation rate for this model is less than about 40%, which is compatible with the maximal value predicted in equation 3.19.
In summary, although introducing new coannihilation partners and s-channel mixed annihilation processes can in principle enhance the dark matter annihilation rate, the concrete possibilities to do so are very limited. For two-partner models, increases of the rate by more than 100% can only be obtained by introducing mediators that have a particularly narrow width, can couple only to different species of coannihilation partners, and annihilate resonantly. Aside from this particularly contrived scenario, it is also possible to simply increase the number of mediator interactions with the dark sector but, as shown in equation 3.19, the resulting enhancement of the dark matter depletion rate will be extremely limited.

t-channel mediators
As described in section 3.1, introducing a new t-channel interaction between X 1 , X 2 , and the SM requires the existence of a mediator M t that is also part of the dark sector, and of two new vertices M t X 1 SM 1 and M t X 2 SM 2 . The existence of these new vertices will not only enable the mixed X 1 X 2 → SM 1 SM 2 interaction but also open new self-annihilation channels for both X 1 and X 2 as shown in figure 4. In order to isolate the effects of the introduction of a new coannihilation partner from the effects related to the mediator itself, we compare the annihilation rates in models with X 1 , M t , and X 2 to the rates in models with either X 1 or X 2 , and the mediator M t .  Figure 9. Ratio of the mixed annihilation rate over the self-annihilation rateR for t-channel models as a function of α max for m M = m X . The scalar and fermion X 8 + X 8 + M 3 models (purple lines) do not appear in figure 8 sinceR is always less than one in these scenarios.
Since the DM and the X particles all have equal masses in the ∆ ≈ 0 region that we are studying here, the t-channel models are described by four parameters: m M , m X , α 1 , and α 2 . For perturbative annihilation cross-sections, the mass dependence is entirely contained in the ratio ρ = m M m X . Similarly to the s-channel procedure described in section 3.3, we computeR(ρ, α max ) defined as R(ρ, α max ) ≡ max α 1 ,α 2 ≤αmax R(ρ, α 1 , α 2 ), (3.20) with R given in equation (3.8). We show the regions of the (ρ, α max ) space where the mixed condition (2.17) is saturated (R = 1) in figure 8 for a subset of the models constructed from equations (3.5) and (3.6). These models have been chosen so that they illustrate all the possible types of behavior encountered in t-channels scenarios. As shown in figure 8, while a few models require large couplings α max 5 α s to reachR > 1, in other models, an enhancement of the dark matter annihilation rate is possible even for moderate couplings. For all the scenarios studied here,R is maximal for m M = m X and the mediator cannot be lighter than m X in order for the dark matter to be stable. Interestingly, for m M ∼ m X the mediator also becomes a coannihilation partner. For the simplified models considered here, however, this particular scenario requires a sizable fine-tuning 1 of the masses of the dark sector particles and we will thus not consider the effects of the mediator coannihilation in the rest of this section.
We now investigate how the ratioR changes as α max increases to arbitrarily high values. In figure 9 we show the ratioR as a function of α max for m M = m X , the value of the mediator mass for whichR is maximal at fixed couplings. We first observe that all tchannel coannihilation models asymptote to a maximumR when α max becomes large. This is rooted in the fact that both the self-annihilation rate induced by the new mediator in the one-partner model and the mixed annihilation rate in the two-partner model scale as α 4 max so the coupling dependence vanishes at large α max . The speed at which this asymptotic value is reached is highly model-dependent. In particular, when X 1 and X 2 have different quantum numbers, as for the X F3 -X F8 -M C3 model, a sharp change in the slope ofR can sometimes be observed for a given m M . This feature indicates a change in the one-partner model giving the largest X self-annihilation rate, and hence in the denominator ofR as shown in equation (3.8). Another important observation is that, for the X C3 + X C8 + M F3 model,R can increase to particularly large values. This behavior is due to the fact that the self-annihilation rates in the one-partner models are p-wave suppressed while the mixed annihilation rate has a s-wave component. Hence,R asymptotes only at values of order 50 for this model. Even for this scenario, however,R remains smaller than two for α max < 4 α s . Moreover, even for large couplings, although the dark matter annihilation rates for the twopartner model will be much larger than for the p-wave suppressed one-partner models, they will remain comparable to the ones obtained in other one-partner models where s-wave selfannihilation is allowed. Hence, we do not expect this particular scenario to challenge the 10 TeV bound found in [6].
We conclude from this analysis that, although models with t-channel mixed annihilation can lead to increased dark matter annihilation rates compared to single-partner coannihilating models, this increase is extremely limited. As in the case of mixed annihilation with a SM mediator explored in section 3.2, obtaining a sizable enhancement of this rate requires either introducing extremely large new physics couplings, or a large number of coannihilation partners. Both scenarios will lead to a rich phenomenology, with a huge diversity of colored particles and interactions that can be investigated at colliders. The next section discusses the expected collider signatures for the three classes of models studied throughout this paper.

Collider phenomenology
In a previous work [6] we studied the possible LHC signatures for simplified dark matter models with a single coannihilation partner. Neglecting the interaction between the dark matter and X, these models can be probed by two types of collider searches: multijet plus missing E T searches targeting the production of a prompt X in association with initial state radiation (ISR) [49][50][51][52][53][54][55][56], and, for scenarios where X and the dark matter are very close in mass, searches for long-lived strongly interacting particles [57][58][59][60][61][62]. Combining the associated LHC constraints with the relic density requirements, the multijet plus missing E T searches are expected to probe the m DM 1 TeV and ∆ 10% region. Searches for long-lived particles, on the other hand, typically probe the ∆ 1% region, that corresponds to multi-TeV dark matter masses [6]. Exploring the remaining region -multi-TeV dark sector particles with moderately compressed spectra -requires both a more powerful collider and new types of detectors, specially designed to look for compressed topologies.
In what follows, we describe how a future 100 TeV collider would constrain the parameter space of single and multi-partner coannihilating models.
While the study detailed in [6] focused on simplified models including only one coannihilation partner, we consider here models with multiple X i . Adding coannihilation partners to a model can have multiple effects on its collider phenomenology. First, the new X i particles will be associated with their own set of constraints from both jets + / E T and long-lived searches. Additionally, the introduction of new particles and vertices is going to modify the pair production cross-section of the X i , possibly introducing new physics X j -mediated processes in addition to the usual QCD processes. Finally, the new vertices associated with the mixed annihilation process will give rise to new production channels for the dark sector particles at colliders in addition to the X i X i pair production. In particular, the mixed annihilation process X i X j → SM SM can be reversed to lead to X i X j production. In models where this process occurs through a new physics mediator, it is also possible to look for the signatures associated to this mediator. In what follows, we evaluate the influence of these different effects for a model with a SM singlet Majorana fermion dark matter DM M and two coannihilation partners X C3 and X M8 . We study the scenario where the mixed X C3 X M8 annihilation occurs through an s-channel quark, keeping in mind that introducing a new physics mediator instead would lead to a richer phenomenology and hence stronger constraints. The dark sector particles in our model are similar to the squark, the gluino, and the bino in SUSY [63,64]. Here, however, we consider only one flavor of squark and model the interaction between the X i and the dark matter using an effective operator.
The parameter space of our two-partner model is spun by the dark matter mass m DM , the relative mass splittings ∆ 1 and ∆ 2 of the two coannihilation partners, the X C3 X M8 q coupling α NP , and the parameters describing the X i DM interaction. Since, as shown in [6], the collider bounds on prompt X i pair-production from multijet + / E T searches have a very weak dependence on ∆ i , we set ∆ 1 = ∆ 2 = ∆ for the rest of this study. As in [6], we treat the X i DM → SM SM coannihilation as a subdominant process for the determination of the dark matter relic density and we model it using an effective operator, with suppression scale Λ. The value of this scale, as well as the structure of the effective interaction will not affect the collider phenomenology of models where X i is prompt; they will, however, affect the lifetime of X i , and hence the constraints from the long-lived particle searches. In what follows, since we will study the phenomenology of our models at a 100 TeV collider, we set Λ = 20 TeV. We use the effective operators introduced in [6] for the X C3 DM SM SM and X M8 DM SM SM interactions, choosing the SM particles to be either quarks or gluons. Note that, with this requirement, since the dark matter is a fermion, one of the products of the X C3 DM coannihilation process will necessary be a gluon. The X C3 DM SM SM effective operator will thus be at least of dimension six and one-loop suppressed. The decay width of X C3 should therefore be extremely suppressed, which will lead to particularly strong bounds from the long-lived particle searches.
We first evaluate the impact of the new X C3 X M8 q vertex described in equation (3.2) on the production rate of X C3 and X M8 at a 100 TeV collider. To this end, we compute the X C3 X C3 , X M8 X M8 , and X C3 X M8 production cross-sections in association with one ISR jet for different values of α NP . We compute these cross-sections at leading order, imposing a mild p T cut of 100 GeV on the ISR jet. We checked that increasing this cut does not affect our conclusions. We inform our choice of α NP by considering the ratio R of the dark matter effective annihilation cross-section in the two-partner model over its value in a one-partner model for ∆ = 0, as defined in equation (3.8). In particular, we choose α NP = 0, 7.0 α s , and 12.9 α s , that correspond to R ≈ 0.57, 1, and 2 respectively, obtained by following the procedure discussed in section 3.2. Note that R min ≈ 0.57 is the minimal possible value for R. The production rates of the coannihilation partners for these different couplings are shown in figure 10 as a function of the mass of the X i . While introducing the new physics coupling always increases the pair-production rate of the X i , this increase remains negligible for the pair-production rate of the X M8 "gluino-like" coannihilation partner. While the "squark" X C3 pair-production cross-section can be increased by up to two orders of magnitude by the introduction of the new couplings, it remains usually lower or similar to the X M8 pair-production cross-section. Finally, the "mixed" X C3 X M8 production rate is comparable to the X M8 pair-production rate for all values of m X . For the multijet plus / E T searches and the compressed topologies that we are studying, all three processes are expected to lead to extremely similar signatures. It is therefore possible to estimate the associated constraints by considering a classic gluino-neutralino simplified model where the gluino production cross-section is multiplied by a factor of a few compared to SUSY. This enhancement, however, would only correspond to an increase of at most a TeV in the mass reach of the multijet search. Hence, the multijet plus / E T constraints on our two-partner model will be similar to the ones associated to the gluino-neutralino simplified model, that have already been computed in the literature [65].
We use the result discussed above to compute the bounds on m X from multijet plus / E T searches for a 100 TeV center of mass energy and 3 ab −1 luminosity using the results derived in [65] for a gluino-neutralino simplified model. In addition, we estimate the bounds associated to future long-lived particle searches using the procedure described in [6], also assuming a luminosity of 3 ab −1 . These different constraints are shown along with the regions of parameter space allowed by the relic density requirement in figure 11. One striking result for the X C3 X M8 model is that the whole R ≤ 2 parameter space allowed by the relic density constraints would be within the reach of long-lived particle searches since the decay rate of X C3 is particularly suppressed. Constraints from these searches can be alleviated if ∆ C3 is increased with respect to ∆ M8 but most of the parameter space of our model is still likely to remain excluded, as such an increase will also considerably tighten the relic density constraints on m DM and ∆ M8 . Note, however, that this result is valid only if X C3 decays via an effective operator. In a SUSY-like scenario where X C3 can decay directly to a quark and a neutralino, for example, constraints from long-lived particle searches will be considerably alleviated. The constraint from multijet plus / E T  Figure 10. Production cross-sections for X C3 X C3 (orange), X M8 X M8 (red), and X C3 X M8 (blue). For each process, we take the new physics coupling to be α NP = 0 (solid), 7.0 α s (dashed), and 12.9 α s (dash-dotted).
searches on the other hand is particularly robust since it does not depend on the structure of the X i DM interaction. For R = 0.57, this search will be extremely efficient and can probe the parameter space down to extremely compressed regions, with ∆ 1%. In these regions, it is expected that the X i will have suppressed decay rates independently from the structure of the X i DM SM SM operator. For R = 1 and 2, the direct searches will probe the parameter space down to ∆ ≈ 2% and 4% respectively. Although a large region remains unexplored, we have to remember that the multijet bounds shown in figure 11 have been computed assuming that the design of the detectors for the future 100 TeV collider would be similar to the one of ATLAS and CMS. In fact, our result shows that designing detectors targeting compressed regions of parameter space would be a crucial step in the search for thermal dark matter [66][67][68][69][70].
Finally, we would like to comment on other smoking gun signatures that could be observed in models with multiple coannihilation partners. First, as noted throughout this paper, in most regions of the parameter space, increasing the number of coannihilation partners in a given model leads to looser relic density bounds only for couplings larger than a few times α s . Such scenarios should therefore be associated with a strongly-coupled sector at energy scales close to m DM . This strongly-coupled sector will typically be associated with a large number of new particles, notably composite states and light Goldstone bosons, that could also be observed at colliders [71][72][73][74]. Additionally, as discussed in section 3, the mixed X i X j annihilation process can be mediated by a new physics particle. This mediator does not need to be close in mass to the dark matter and will therefore give rise to characteristic collider signatures. In particular, for s-channel models, the mediator can be light and decays to two SM particles. It could therefore be a privileged target for resonance searches [75][76][77][78][79][80][81][82][83][84][85][86][87], with signatures similar to the ones studied in [7,8]. In the D MM Co llid er Jet s Li  Figure 11. Example model for colored coannihilation which shows the general interplay among collider limits and relic abundance constraints. We compare two models with either a single X C3 (blue) or a single X M8 (orange) coannihilation partner to a model with both X C3 and X M8 at the same time. We show three scenarios for the two partner model: one where the mixed annihilation is absent (green) and two other scenarios with R = 1 (red) and R = 2 (purple). Imposed are the exclusion limits from a prospected 100 TeV hadron collider for the single partner models.
t-channel, the decay mode of the mediator will be similar to the one of the X i but since the mediator does not have to be close in mass to the dark matter, multijet plus / E T searches should now be much more powerful.

Conclusions
Confronting experimental limits to relic density constraints has been a major tenet of the search for thermal dark matter during the last few years. In particular, for coannihilating dark matter scenarios, using the simplified model formalism has allowed to translate results from the LHC searches into limits on the mass of the dark matter [88][89][90][91][92][93][94][95][96] and on its splitting with its coannihilation partner X. These results in turn have been used to motivate both building a 100 TeV collider and develop new detector technologies, targeted towards models with particularly compressed spectra. In this study, we evaluated the robustness of these simplified model results against introducing new coannihilation partners.
Throughout this paper we focused on models with colored coannihilation partners that are charged under QCD although our results can be straightforwardly extended to other models with unbroken gauge groups. We find that, in general, increasing the number of coannihilation partners in these models lead to tighter relic density constraints, and hence to lower dark matter masses, in most regions of the parameter space. The few regions in which adding new coannihilation partners allows for larger dark matter masses correspond to scenarios where the dark matter effective annihilation rate is dominated by "mixed" annihilation processes between different species of coannihilation partners. We showed that these different regions can all be characterized by focusing on models with only two coannihilation partners.
For a representative subset of these models, we considered three possible classes of mixed annihilation processes: mediated via an SM quark, and s-and t-channel via a new physics mediator. In the first class, the mixed annihilation rates dominate only when the new physics coupling of the SM to the coannihilation partners is much larger than the strong coupling. In the second class, introducing a new coannihilation partner can lead to a significant increase of the dark matter annihilation rate near the resonant annihilation region, where the mediator mass is close to the sum of the masses of the coannihilation partners. This increase, however, only takes place when the mediator of the mixed annihilation process has a particularly narrow width. Finally, mixed annihilation processes happening in the t-channel can dominate over similar self-annihilation processes only if the latter are associated to particularly small rates, either being p-wave suppressed or involving color triplet particles with a small spin. Since such suppressed processes lead to relatively tight relic density bounds, of around a TeV, on the dark matter mass raising the dark matter annihilation rate for these models is unlikely to affect the current estimates of the upper bound on the energy scale of the coannihilating models, which lies around 10 TeV. This upper bound can therefore be loosened only in very specific models involving resonant s-channel annihilation processes, or in models involving particularly large new physics couplings, that are likely to exhibit a non-perturbative dynamics.
We finally studied the possible signatures associated with models with multiple coannihilation partners at a future 100 TeV collider. We focused on a simple SUSY-like scenario where a scalar color triplet and a Majorana color octet coannihilate with a Majorana fermion dark matter candidate. We showed that, for this scenario, the interaction between the triplet and the octet only leads to order one changes in the total pair-production rate of these new particles. We also showed that the different coannihilation particles in a given model can be associated to complementary bounds -with one particle having a large pair-production rate and the other one being long-lived for example -that would lead to a spectacular reach. Finally, we would like to emphasize that, even in regions where the dark matter annihilation rate is increased compared to the one-partner models, a 100 TeV collider could probe mass splittings between the dark matter and its coannihilation partners down to a few percent. In order to explore the remaining region, it is essential to develop new detectors, experiments and analysis techniques, specially geared towards these small mass splittings.
We showed in this paper that, even in non-minimal scenarios, dark matter models with an arbitrary number of coannihilating partners generally cannot have dark matter particles heavier than about 10 TeV. Exceeding this energy scale requires introducing either new particles and vertices or particularly large couplings, which would imply a particularly rich and diverse phenomenology at colliders. Exploring this phenomenology should therefore be a crucial element of the dark matter search program at a future 100 TeV machine. Finally, we would like to emphasize that the methodology detailed throughout this paper can be straightforwardly adapted to study other extensions of the current simplified dark matter models, such as scenarios with multiple dark matter candidates or new unbroken gauge groups. In a time where minimal scenarios are being increasingly cornered, our approach then provides an effective toolkit to comprehend complete models and thus guide the design of the next generations of colliders to allow them to say the final word on thermal dark matter.

A Condition on the effective annihilation cross-section
In this appendix we consider a general model with one dark matter candidate and N coannihilation partners as described in section 2.1.3. Recall that the corresponding effective annihilation cross-section as defined in equation (2.13) is where we now wrote the double sum explicitly. Here, we show that requiring all submodels of the form DM + X i + X j to satisfy (2.12) is a sufficient condition for σ eff X 1 ···X N to be smaller than max σ eff X 1 , . . . , σ eff X N .
Without loss of generality, we assume that σ eff We then saturate the constraint in equation (2.12) for each σ mix X i X j to obtain (A.2) Denoting by σ eff max the maximal effective annihilation cross-section for models with only one of the X i we can then write In order for the effective annihilation cross-section for the complete model to be maximal σ eff max , we need the ratio to be at most one. We note by explicit computation that which implies that R has a stationary point at R = 1 that can be approached from either R > 1 or R < 1 but not both simultaneously. Now if we can show that one point in the g X i space has R < 1 we then know that for all points R ≤ 1. In fact it can be shown that for the point g X i = g DM we always have R < 1 for N > 1. Hence, when all the DM + X i + X j models verify condition 2.12, the dark matter effective annihilation rate is also decreased when all the coannihilation partners are present simultaneously.

B Sommerfeld corrections
In our analysis, we take the Sommerfeld corrections into account when computing the annihilation cross-sections of the colored coannihilation partners. We base ourselves on [97] and on a previous work [39] which focuses on the self-annihilation rates of colored particles into gluons and quarks through the strong interaction. Here we extend this work to include the Sommerfeld corrections for the mixed annihilation rates and the new physics contributions to the self-annihilation rates stemming from the interactions shown in the Lagrangians (3.2), (3.4), (3.5) and (3.6) and depicted in figures 3 and 4. We adopt here the same strategy as in [39], decomposing the product of the SU (3) representations of the initial state particles into a sum of color eigenstates, each of them associated to a Coulomb potential. This procedure allows to derive particularly simple expressions for the Sommerfeld corrections to the the self-annihilation rates of colored particles using CP conservation and the symmetry properties of their color representations. For mixed annihilation rates, which involve initial state particles of different types, the situation is more complicated and we need to consider each process individually. In the rest of this appendix we use the formalism and notation of [39]. Moreover, we update the Mathematice package [48] to compute Sommerfeld corrections for the processes discussed here.
Although non-perturbative effects like bound state formation [12] can occur alongside Sommerfeld corrections, we focus solely on the latter for this study. As shown in [6,12], the effect of bound state formation on the effective annihilation cross-section is generally milder than the Sommerfeld corrections. When taking ratios of cross-sections as for the mixed condition derived in section 2.1 these non-perturbative effects partially factor out between the one and two-partner models. We will therefore not consider bound state formation throughout this work.
The sizes of the non-perturbative corrections to the mixed annihilation rates heavily depend on the quantum numbers of the coannihilation partners involved as well as the structure of the annihilation process. In this section we adopt the strategy detailed in [39], computing the Sommerfeld corrections separately for each product of two diagrams that enters the total cross-section, that is, squared and interference terms. As discussed in section 3 mixed annihilation can proceed through either a SM quark, an s-channel mediator or a t-channel mediator. Each scenario involves a different set of diagrams, each of them associated to a specific color structure.
We first notice that there is a wide range of possibile configurations for the mixed annihilation diagrams since the annihilating particles can transform as triplets, sextets, or octets under the strong gauge group. We thus need to calculate the QCD decomposition of all the possible pairs of colored initial states that can annihilate into colored SM particles. All the relevant initial state color combinations for either self-annihilation or mixed annihilation as well as the conjugate expressions. Analogous to [9,39] the QCD potential for an initial state whose particles have representations R and R is decomposed into channels of definite color as The coefficient α Q can be calculated for all the decompositions in equation (B.1) and is given in the following table. In this table a blank entry indicates that the color decomposition of the representation product R ⊗ R does not contain Q. We note that the values of α Q are independent of the considered annihilation processes and their diagrammatic structures.
The remaining necessary step to calculate the Sommerfeld corrections is to compute the size of the contributions of the different Q initial states to the total annihilation cross-section. As mentioned earlier in this section we need to treat each product of two diagrams entering in the cross-section separately. In section 3.1 we detailed which diagrams are involved in the mixed annihilation processes. The vertices involved in these diagrams are associated with the following color factors where the indices give the color representations of the particles involved. Since Sommerfeld corrections for self-annihilation cross-sections have been computed in [39], we will focus here on the mixed annihilation rates. However, in models involving a SM mediator or a t-channel NP mediator, allowing for mixed annihilation also opens new channels for the self-annihilation of the X i into SM quarks. The corresponding diagrams interfere with the existing QCD s-channel annihilation process and we need to compute the size of the Sommerfeld corrections for the associated new contributions.
For each annihilation process, we write the total cross-section as the sum over all the color This allows us to write the Sommerfeld corrected cross-section in the notation of [39] as Note that if α Q is negative we have Sommerfeld enhancement and if it is positive we have a Sommerfeld reduction of the annihilation rate in the specific color channel Q.
We first look consider processes occuring via a s-channel NP mediator. These cases are particularly simple to treat as the color representation of the initial state has to be the same as the one of the mediator, and we therefore obtain β Q = 1 for a mediator with charge Q. Then the single Sommerfeld correction factor α Q can be read of table (B.3). The s-channel NP models do not induce self-annihilation processes for the coannihilation partners so we do not need to compute any additional corrections for these models. A similar simplification occurs for the interference of the new physics self-annihilation processes in SM and tchannel models with the QCD self-annihilation into SM quarks. Since the latter proceeds through an s-channel gluon these terms are always associated with β 8 = 1, and non-octet initial states do not contribute. Finally, the SM mediated mixed annihilation of X 1 and X 2 into a SM quark and a gluon involves one s-channel and two t-channel diagrams. For the s-channel amplitude squared and any of its interferences with the other diagrams the previous argument applies and we have β 3,3 = 1 depending on whether the mediator is a quark or an antiquark.
The only remaining contributions to the annihilation cross-sections are the terms involving only products of t and u-channel amplitudes, that appear in the SM and t-channel NP models. The corresponding processes can be uniquely represented by R 1 R 2 → S 1 S 2 where R i and S i are the color representations of the initial and final states respectively and we denote the color of the mediator with M. The color flow for these types of diagrams is defined in figure 12. We denote the color structure of the diagram as C R 1 R 2 ,M,S 1 S 2 , keeping in mind that |C R 1 R 2 ,M,S 1 S 2 | 2 = |C R 2 R 1 ,M,S 2 S 1 | 2 . For the initial state representations and their color decompositions described in equation (B.1) there are many possible combinations for C. In the next tables we list the cross-section decomposition coefficients β Q for the models introduced in section 3. For the SM-mediated models discussed in section 3.1.1, the squared and mixed products of the three leftmost diagrams in figure 3 are associated with the following coefficients In these tables some of the representations Q appearing in the decomposition (B.1) are not present since they do not appear in the color decompositions of the SM final states and therefore cannot contribute to the annihilation processes. Since the diagram in figure 12 also has u-channel variants we use the superscripts t and u to indicate the respective channel whenever the color flow of the interference term is ambiguous. Now, using equation (B.6) and by reading of the values for α Q in table (B.3) and for β Q in tables (B.7) and (B.8), we can calculate the Sommerfeld corrections for all processes considered in this work.