LHC Phenomenology of Dark Matter with a Color-Octet Partner

Colored dark sectors where the dark matter particle is accompanied by colored partners have recently attracted theoretical and phenomenological interest. We explore the possibility that the dark sector consists of the dark matter particle and a color-octet partner, where the interaction with the Standard Model is governed by an effective operator involving gluons. The resulting interactions resemble the color analogues of electric and magnetic dipole moments. Although many phenomenological features of this kind of model only depend on the group representation of the partner under SU(3)$_c$, we point out that interesting collider signatures such as $R$-hadrons are indeed controlled by the interaction operator between the dark and visible sector. We perform a study of the current constraints and future reach of LHC searches, where the complementarity between different possible signals is highlighted and exploited.

the SU(3) c gauge group [25]. However, as we will discuss below, important features for collider phenomenology are indeed reliant on the interaction between the dark sector and the SM.
In this paper, we consider the SM augmented by a dark sector constisting of a DM particle and a nearly-degenerate colored state, in the adjoint representation of SU(3) c . This dark sector communicates with the SM via a dimension-5 effective operator (the validity of effective theories for DM searches has been widely discussed in the literature, see e.g. Refs. [29][30][31][32][33][34][35][36][37]). Such a scenario is particularly interesting because the colored partner could hadronize in bound states like ordinary quarks and gluons. In a supersymmetric context this is a well-known possibility, and such bound states, originally introduced in Ref. [38], are called R-hadrons. We use here the same terminology, although our considerations do not assume any underlying supersymmetry. For more recent papers about R-hadrons, see e.g. Refs. [39][40][41][42][43]. In addition, since the decay of the colored partner is governed by a suppressed non-renormalizable operator, such a bound state can easily travel macroscopic distances and leave tracks in the collider detector.
The paper is organized as follows: in section 2, we introduce the model and discuss some of its features and implications; in section 3, we consider LHC constraints derived from monojet and R-hadron searches, focusing on the interplay between them. Finally, we conclude in section 4.

Model
Dark matter, despite being neutral, can be coupled to colored Standard Model particles. In order to allow such a coupling, colored particles within the dark sector are required [28,[44][45][46][47][48][49][50][51]. In this work, we consider an extension of the minimal scenario, where the DM particle χ 1 is accompanied by a slightly heavier partner χ 2 . We denote the masses of these particles by m 1 and m 2 ≡ m 1 +∆m, respectively. Both χ 1 and χ 2 are Majorana fermions.
At the renormalizable level, scalar or fermionic partners in the fundamental representation of SU(3) can be responsible for the coupling of DM with the SM quarks [28,52,53]. If, instead, we are to consider a coupling to gluons, the lowest dimensional operator has D = 5 and involves a colored partner χ 2 in the adjoint representation of SU (3). If we denote the dark matter particle by χ 1 , the free Lagrangian for the dark sector is: with a being the color index in the adjoint representation. The coupling to gluons can be attained via effective operators mimicking the (chromo-) electric and (chromo-)magnetic dipole moments, as follows: where σ µν = i/2[γ µ , γ ν ] and G a µν is canonically normalized. The two operators in Eq. (2.2) give rise to similar phenomenology and no interference effect arises in any of the observables we study in this paper. Therefore, for simplicity, we limit ourselves to study only the operator with d χ in the rest of the paper.
Effective operators describing dipole moments typically arise after integrating out heavy particles of the underlying ultraviolet theory at loop-level. If this is the case, the operator in Eq. (2.2) should be further suppressed by α s /(4π), and so the importance of higher-order operators may not be negligible. A more complete theoretical analysis of the origin of the interaction in Eq. (2.2) and of the role of higher-order operators is left to a future work. ij SM SM qq g g Table 1. Different contributions to the effective cross-section σv χ i χ j →SM SM . The QCD coupling is denoted by gs, while v is the relative velocity in the χiχj center-of-mass frame.
The interactions of the dark sector with the SM particles are then described by the parameters {m 1 , ∆m, d χ }. In particular, we require that d χ 1: this interaction term, in fact, could be written as an effective term suppressed by 1/Λ, with Λ being the scale of some underlying new physics. It is then natural to formally identify d χ ∼ m 1 /Λ, which has to be small in order for the effective theory to be reliable. For the values of d χ considered in our analysis, the energy scales of the processes of interest are always well below the operator scale Λ, thus ensuring we are in the regime of valid effective field theory.
The simplest process leading to the decay of χ 2 is χ 2 → χ 1 g, whose width, at leading order in ∆m/m 1 , is: Since d χ is required to be small, we would naturally expect χ 2 to be a long-lived particle with lifetime on the detector timescale, as will be explored later.

Relic density
A first constraint on the parameter space can be obtained by requiring that the model reproduces the observed dark matter abundance Ωh 2 = 0.1194 [1]. Such a relic density is determined by processes of the form σ(χ i χ j → SM SM). The expressions for the corresponding cross sections, at leading order in ∆m/m 1 and m f /m 1 (where m f is the mass of a generic SM fermion f ), are shown in table 1, although the complete expressions, which can be found in Appendix B, have been used in the calculations. In order to determine the relic density predicted by this model, two modifications to the standard procedure have, in principle, to be taken into account: first, if the mass splitting between χ 1 and χ 2 is small compared to their masses, co-annihilations must be included [26,54]; second, due to the color charge of χ 2 , Sommerfeld enhancement (introduced below) modifies the value of σ(χ 2 χ 2 → SM SM). As far as the co-annihilations are concerned, the effective cross-section which determines the observed abundance of DM in the universe is: where α ≡ g 2 /g 1 (1 + ∆m/m 1 ) 3/2 e −x∆m/m1 , x ≡ m 1 /T , σv ij ≡ σv χiχj →SM SM and g i is the number of degrees of freedom of χ i . The relic abundance is then related to this effective cross-section as: where g * is the number of relativistic degrees of freedom at the freeze-out temperature T F , determined by the implicit equation: with n F = g 1 (1 + α) being the effective number of degrees of freedom of the system (χ 1 , χ 2 ). In the following, we take g * = 106.75 as a reference. As already mentioned, Sommerfeld enhancement also plays an important role in the determination of the relic abundance [55][56][57]: Sommerfeld enhancement is a non-perturbative effect due to the exchange of soft gluons betweeen the colored particles in the initial state. This is therefore relevant for the self-annihilation of χ 2 . Model independent discussion of this effect can be found in Refs. [25,27]. These analyses assume that the relic DM density is dominated by QCD, remaining agnostic about the particular phenomenology deriving from the new BSM coupling. This is a reasonable assumption for the model we consider since, as already stated, it is natural (and indeed necessary) to assume that d χ 1. Co-annihilations where the DM annihilation cross-section contributes negligibly to the relic density have been recently analyzed in Ref. [58] in the more general context of sterile co-annihilations.
When the final state is characterized by a single representation Q, the Sommerfeld-corrected cross-section is σ Somm = S (C Q α s /β) σ Pert , where S is the non-perturbative correction depending on the final representation (through the Casimir element C Q ) and on the velocity of the particles β. If, on the other hand, we have more than one possible final state representation, we need to consider the decomposition R ⊗ R = ⊕ Q Q, where R and R are the initial state representations (in our case R = R = 8) and Q's are the final state ones. Each representation Q gives a contribution to the total cross-section and has its own value of C Q . After group decomposition, the final result is given by eqs. (2.24, 2.25) of Ref. [25].
As a result, the contour yielding the correct relic density, with and without the inclusion of such a non-perturbative effect, can be found in fig. 1. In this plot, we only consider the dominant contributions from QCD self-annihilations, not including the sub-leading contributions of processes proportional to d χ , which will be negligible if d χ 1. This is actually well motivated from the previous discussion about the magnitude of d χ .
While the relic density is dominated by QCD processes, we see from eq. (2.3) that the decay length instead depends quadratically on d χ . Therefore the smallness of d χ leads to macroscopic decay lengths, which are an interesting feature we will use in our analysis. From here on, we fix the mass splitting ∆m as a function of the mass of the DM candidate m 1 , using the Sommerfeld corrected curve in fig. 1. This imposes the correct relic density for all points in parameter space that we consider. As a consequence, the decay length now only depends on the mass of χ 1 and the coupling d χ . The full numerical results for the decay length can be found in fig. 2, where we show contours of the proper decay length of χ 2 .

Departure from chemical equilibrium
The co-annihilation paradigm, described in section 2.2, implicitly assumes that chemical equilibrium is maintained until the DM freeze-out. Under particular conditions, however, it is possible that this assumption might not be valid [28]: this can happen, for instance, when the relic abundance is dominated by a SM (here, QCD) coupling: in this case, the coupling characterizing the BSM physics (d χ in the case at hand) remains unconstrained. This translates into the fact that very small values for such a coupling are in principle allowed, leading to a possible breakdown of chemical equilibrium.
The important ratios to evaluate are Γ χiχj /H, where Γ χiχj generically represents the rate of a process involving χ i and χ j : it can be the scattering χ 2 χ 2 → χ 1 χ 1 , the decay χ 2 → χ 1 g, the conversions χ 2 g → χ 1 g and χ 2 q → χ 1 q, as well as all the inverse reactions.
The rates of decay and conversion are proportional to d 2 χ , while that for the scattering is proportional to d 4 χ ; therefore this latter process is expected to be the most sensitive to d χ , and is therefore expected to have the smallest rate. When the largest of these rates Γ χiχj /H 1, the assumption of chemical equilibrium (which eq. (2.4) relies on) ceases to be valid. If this is the case, a numerical integration of the complete set of Boltzmann equations, including conversions, is necessary. The ratios Γ/H for these three processes (in the direction χ 2 → χ 1 ) are shown in fig. 3.
Since the rate corresponding to σv χ2 g↔χ1 g ∝ g 2 s d 2 χ turns out to be the dominant contribution, scatterings with gluons are ultimately responsible for maintaining chemical equilibrium.
In order to test the possible breakdown of chemical equilibrium before and during freeze-out (where eq. (2.4) has its validity), the ratio Γ χ2 χ1 /H ≡ H −1 n σv χ2 g↔χ1 g can be investigated in the region 20 x F 30. From fig. 3, we see that in this region Γ χ2 χ1 /H ∼ 10 4 for d χ = 10 −6 . We can therefore estimate the breakdown of chemical equilibrium to occur when: This simple scaling argument is actually in agreement with the explicit result shown in fig. 3. In the following, we therefore assume d χ 10 −8 , in order to be in the regime of chemical equilibrium.

LHC searches
In this section we analyze the constraints on the model coming from the two most important channels: R-hadrons and monojet. In principle, it would also be possible to have limits from dijet-resonance bounds coming from the production and fragmentation of a bound state of two χ 2 particles, similar to a gluinonium. Since this results in rather weak constraints, we have described it in Appendix A.

R-hadron constraints
The color charge of the χ 2 particle implies that it can hadronize with SM particles on the detector timescale, forming particles analogous to the R-hadrons in supersymmetry. If stable on a detector timescale, these colorless composite states can be detected via an ionization signature as they travel through the detector at speeds significantly less than the speed of light.
We apply ATLAS constraints on the χ 2 production cross-section from Ref. [39], which searches for R-hadrons at √ s = 13 TeV with 3.2 fb −1 of data. The relevant constraints are those on gluinos, since χ 2 is a color octet. We also consider an approximate high-luminosity (HL) projection of these limits to L = 3000 fb −1 , using the procedure outlined in Ref. [24], applied to the ATLAS analysis. The relevant results are the background counts in Table 3 of Ref. [39] for the gluino search, which we rescale with the increased luminosity. We assume that the same efficiencies of Table 3 apply to the HL bounds. It should be noted that in the HL regime the results are limited by systematics rather than statistics. The signal simulations are the same for the current luminosity and for higher luminosities.
In order to simulate the pair production of χ 2 particles at parton level we have used Mad-graph5 aMC@NLO [59], where the model has been implemented using FeynRules [60], and apply the R-hadronisation routine from Pythia 8.230 [61]. The probability of each χ 2 being stable at least up to the edge of the ATLAS calorimeter is given by where calo = 3.6 m is the transverse distance to the edge of the calorimeter and we defined T = p T 2 /(m 2 Γ χ2 ). This probability is applied on an event-by-event basis to find the effective crosssection of events yielding at least one R-hadron. This relies on the assumption that the lifetime of the resultant R-hadron is at least as long as the unhadronized χ 2 lifetime. Following Ref. [39], we assume that 90% of the χ 2 form charged R-hadrons.
Contours showing the relationship between this effective production cross-section, m 2 and d χ are shown in Fig. 4, along with current and projected future ATLAS limits on the cross-section.

Monojet
A generic particle physics model for DM is usually sensitive to so-called 'monojet' searches, where DM produced in a collider recoils from a high-energy jet, leaving a large missing energy (E miss T ) signature as it passes through the detector without interacting [9,[62][63][64].
For the chromo-electric model the production processes leading to the monojet signature are of the form p p → χ i χ l j with i, l ∈ {1, 2}. Since d χ 1, the leading contribution will be from the QCD-mediated production channel p p → χ 2 χ 2 j, since all other terms are proportional to powers of d χ . In this regime, the relic density profile shown in fig. 1 will apply. We apply the latest monojet constraints from ATLAS [9], which searches for events with large missing energy and at least one high-energy jet, with center-of-mass energy of 13 TeV and integrated luminosity of 36.1 fb −1 .
Events are required to satisfy the conditions E miss T > 250 GeV, leading-p T > 250 GeV and also |η| leading−jet < 2.4. In addition, a maximum of four jets with p T > 30 GeV and |η| < 2.8 are allowed, and the condition ∆φ(jet, p miss T ) > 0.4 must be satisfied for each selected jet. The analysis then uses ten different signal regions, which differ from each other by the choice of cut on E miss T : in particular, the weakest one is denoted (for the inclusive analysis) by IM1, and requires E miss T > 250 GeV; while IM10 requires E miss T > 1000 GeV. We simulate events at parton level using Madgraph5 aMC@NLO [59], then apply the same cuts as Ref. [9]. In models where a colored partner is produced at the LHC, monojet constraints will only apply if the colored partner decays promptly, i.e. within the beamline radius, beam = 2.5 cm. Otherwise, if it enters the detector material, it will form an R-hadron within a very short timescale, roughly Λ −1 QCD ∼ 10 −24 s. We take this into account by considering the probability for each particle χ 2 to decay with transverse decay length T less than d beam [27]: where T = p T 2 (i)/(m 2 Γ χ2 ) is the transverse distance traveled by χ 2 in an event i. Each event is weighted by this probability in order to find the effective cross-section where χ 2 decays promptly, before forming an R-hadron. We assume here that all the colored particles reaching the detecting stage hadronize.
In order to obtain a limit on the number N NP of new physics events, for both current and future luminosities, we apply a χ 2 analysis with 95% CL with unit efficiency and acceptance, according to [25]: where the error on the SM background is assumed to be normally distributed.
To find the strongest constraint, we consider the different signal regions from Ref. [9], differing by the cut on E miss T . For a given value of m 1 , we use the ratio between our simulated cross section and the bound from the ATLAS paper and find that the strongest bound comes from IM9 (which requires E miss T > 900 GeV) as can be seen in fig. 5. It should be noted that changing the mass varies both the value of the cross section and the kinematic distribution of the particles, so that the results from fig. 5 cannot be trivially recast into a bound on the mass.
For our optimal bin, the number of events in this signal region is: Then the cross section of new physics (NP) has to satisfy the constraint σ N P < 2.3 fb for L = 36.1 fb −1 . Using this value and the procedure outlined earlier in this section, we find a lower bound on the mass of the DM of 860 GeV for d χ 3 × 10 −7 . Full results are shown as the blue lines in Fig. 6. For smaller values of d χ , χ 2 begins to travel into the detector and form R-hadrons before decaying, as discussed earlier in this section.
We extrapolate the monojet bound from Ref. [9] to higher luminosity by considering the statistical and systematic uncertainties separately. The relative statistical error scales with the inverse square root of the number of events (and hence of the luminosity); on the other hand, it is generally Figure 5. Ratio between the cross section from our model and the bound from Ref. [9] in the case of m1 = 860 GeV as a function of the inclusive regions.
not straightforward to predict how the relative systematic uncertainty will evolve with the luminosity. For this reason, we parametrize it in general as δ sys (L 2 ) ≡ r δ sys (L 1 ). Using the published upper bound on the cross-section of new physics, σ N P , at a luminosity L 1 , we can then estimate the corresponding upper bound on σ N P at a different luminosity L 2 as: We carry out this HL projection to 3000 fb −1 , in an optimistic scenario where systematic uncertainties have been cut to half the current values. In this case, we find a limit on m 1 from monojet of 1020 GeV. Results are shown as the dashed blue curve in Fig. 6. For reasons of completeness, we also considered the extreme cases in which the systematics will be unchanged with respect to their current value and the case in which they will be completely negligible, getting respectively the bounds m 1 > 900 GeV and m 1 > 1250 GeV.

Comparison between different searches
The chromo-electric dipole dark matter model has been analyzed in the light of different LHC signals, namely monojet and R-hadrons.
Since the different LHC signals are most effective in different regions of parameter space, it is important to understand the interplay between them. A first noteworthy feature is that the monojet analysis is insensitive to the value of the BSM coupling d χ in most of the parameter space, but at some point this search becomes ineffective due to the fact that the colored partners are forming R-hadrons, rather than decaying to the DM. Since the observed events in this region of parameter space are the R-hadrons, the search for these states becomes the one giving the most stringent bound, as can be seen from fig. 6, where the result of the previous section are summarized. Note that the R-hadron results are shown here in terms of m 1 , since relic density fixes ∆m for given m 1 and d χ .
The complementarity of the searches emerges from the fact that for d χ 3 × 10 −7 the most stringent bound is given by monojet searches, while for d χ 3 × 10 −7 the R-hadron search gives the best result.
Furthermore, the different analyses are affected by different errors, meaning that increased luminosity has a distinct effect on each of them. This suggests that a high-luminosity projection might tell us which of these searches will become more interesting in the future. This as well is shown in fig. 6, where the role of higher luminosity in probing the parameter space is manifest. Figure 6. Current (solid) and foreseen future (dashed) status of the parameter space as excluded by monojet and R-hadrons searches, in blue and orange, respectively. The mass splitting ∆m is fixed, for given m1, by the relic density as shown in fig. 1. The assumptions made regarding the scaling of the error at high luminosity can be found in the text.
As a side remark, we also checked the indirect detection limits by applying bounds on the self-annihilation rate derived from cosmic-antiproton fluxes [65] to our model. An upper bound on d χ , as a function of m 1 , is obtained from the upper bound on the annihilation cross section σ(χ 1 χ 1 → gg). The bounds were found to be very weak compared to those from collider searches, and in a region where the requirement d χ 1 was not satisfied, e.g. the upper limit on d χ was found to be d χ ≤ 0.2 for m 1 = 1 TeV and d χ ≤ 1 for m 1 = 5 TeV.

Conclusions
In this paper, we have explored a remarkable possibility for DM phenomenology at the LHC: the combination of monojet and R-hadron searches. We performed our analysis using a simple effective operator of dark matter, as a case study giving rise to such a situation.
Since the cosmological abundance is dominated by QCD interactions, the coupling of the effective operator d χ is not fixed by the relic density requirement, but it remains a free parameter. The only assumption we make is d χ 1, in order for the effective theory to be reliable. If d χ is small enough, the chemical equilibrium can break down before the dark matter freezeout. We analyzed such a situation and concluded that for the parameter space of interest for LHC searches (d χ 10 −8 ) there is no need to take into account the breakdown of chemical equilibrium.
Our main analysis consisted of the combination of monojet and R-hadrons searches, and found the regions of the (m 1 , d χ ) parameter space excluded by current searches (see Figure 6): while current monojet is able to exclude all points of the parameter space for m 1 900 GeV for d χ 3 × 10 −7 , current R-hadron results are instead able to constrain the parameter space for larger masses, but smaller couplings. This complementarity is maintained in higher-luminosity projections.
These results show once more the importance of finding complementary phenomenological signatures and the power of their combination in strengthening the reach of LHC searches for dark matter.  Figure 7. left: ATLAS excluded (shaded) and projected future limit (dashed) on dijet production crosssection assuming 100% acceptance and branching ratio into dijets, along with theoretical production crosssection (dot-dashed). right: Corresponding current excluded (shaded) and future limit (dashed) on the Br × A. Maximum possible value Br × A = 1 shown as dotted line to guide the eye.

A Estimate of dijet constraints
If a pair of χ 2 are produced near rest, then rather than promptly decaying or forming R-hadrons, they can combine to form a QCD bound state analogous to gluinonium. The χ 2 χ 2 production rate is controlled only by QCD processes, and so should be modelindependent. Therefore we use the production rate of gluinonia as calculated by Ref. [66] as an estimate of the production rate of χ 2 χ 2 bound states. The strongest constraints on this channel come from limits on the dijet resonance production cross-section. We use model-independent constraints from ATLAS [67], taking the conservative choice of assuming a narrow Gaussian width. We also use these limits to estimate the future high-luminosity constraints at 3000 fb −1 using the same method as in previous sections. In order to evaluate the bounds, we have worked under the assumption that the fitting function for the dijet mass distribution used in [67] is still suitable at higher luminosities and we have also considered the case in which the systematic uncertainties are unchanged in the projection.
These constraints are shown in fig. 7 (left), along with the theoretical production rate from Ref. [66]. The limits assume that both the branching ratio to dijets (Br) and the acceptance (A) are 1. In reality Br × A < 1, and the limits are weakened proportionally. In fig. 7 (right), these same results are visualized as an upper limit on Br × A, above which the model is ruled out.
For values of Br×A greater than the maximum theoretical value of 1, the model is not currently constrained by dijets. For m 2 650 GeV, dijet constraints rule out the model only if Br × A is greater than around 0.1 -1.
This constraint is conservative for two reasons: first, the production rate calculation is performed at 14 TeV while the constraints are at 13 TeV. The true 13 TeV production cross-section will be slightly smaller than shown; second, we have taken the dijet constraints assuming a narrow Gaussian width. A broader width weakens the constraints, as seen in Fig. 5 of Ref. [67]. In conclusion, while dijet searches do not strongly constrain the model at the moment, they may be an interesting channel to study with future data.

B Analytical expressions for the differential cross sections
Here we list the analytical expressions for the differential cross sections, separately for each process. The expressions are simplified by using the quantities E 1 = m 2 1 + p 2 and E 2 = m 2 2 + p 2 , which are the energies of the incoming particles χ 1 and χ 2 . The expressions for the cross sections are referred to the center of mass frame (CM), so that p and θ must be interpreted as the momentum of the incoming particles in the CM frame and the angle between incoming and outcoming momenta in the CM frame. The decay to massive quarks is also considered, here m q refers to the mass of the final quarks and the value for the scattering cross section refers to just one flavor.
The self-annihilation of the DM particles into gluons proceeds via t and u channel exchanges of χ 2 .
The annihilation of the partners into gluons proceeds via five different channels: t and u exchange of χ 1 , whose amplitude will be proportional to d 2 χ , s-channel gluon interaction, and t and u exchange of χ 2 , whose amplitude will be proportional to g 2 s . Hence the total cross section will be the sum of three terms proportional to g 4 s , g 2 s d 2 χ and d 4 χ .