A Reappraisal on Dark Matter Coannihilating
with a Top/Bottom Partner
WaiYee Keung ^{*}^{*}*keung@uic.edu, Ian Low ^{†}^{†}†ilow@northwestern.edu, Yue Zhang ^{‡}^{‡}‡yuezhang@northwestern.edu
Physics Department, University of Illinois at Chicago, IL 60607, USA
Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA
High Energy Physics Division, Argonne National Laboratory, Argonne, IL 60439, USA
Abstract: We revisit the calculation of relic density of dark matter particles coannihilating with a top or bottom partner, by properly including the QCD boundstates (onia) effects of the colored partners, as well as the relevant electroweak processes which become important in the low mass region. We carefully set up the complete framework that incorporates the relevant contributions and investigate their effects on the cosmologically preferred mass spectrum, which turn out to be comparable in size to those coming from the Sommerfeld enhancement. We apply the calculation to three scenarios: binostop and binosbottom coannihilations in supersymmetry, and a vector dark matter coannihilating with a fermionic top partner. In addition, we confront our analysis of the relic abundance with recent direct detection experiments and collider searches at the LHC, which have important implications in the binostop and binosbottom scenarios. In particular, in the binostop case recent LHC limits have excluded regions of parameter with a direct detection rate that is above the neutrino floor.
NUHEPTH/1702
1 Introduction
Understanding nature of dark matter (DM) in our universe is of tremendous importance and interest to particle physics and cosmology. If the DM is made of a new particle, a new theory beyond the standard model (SM) will be called for.
This work is devoted to studying a scenario where the DM particle is accompanied in the mass spectrum with a heavier partner which is colored and participates in the strong interaction. Both are odd under a symmetry which is good enough so that the DM is cosmologically stable. Their quantum numbers are chosen such that a renormalizable interaction is allowed among the DM, its partner and a third generation SM quark. We will focus on the class of models where the DM is a gauge singlet under the SM so that if its partner is too heavy and decoupled it has no efficient annihilation channels and will be thermally overproduced. If the partner mass is nearby, the DM could obtain the observed relic abundance in the early universe through the socalled coannihilation mechanism [1]. During the thermal freeze out the DM particles are most efficiently depleted by converting into the partner particles which can annihilate away through strong interactions. By adjusting the DMpartner mass difference (typically less than 100 GeV), the correct relic abundance can be explained for a wide range of DM mass from hundreds of GeV to multiple TeV scale. The coannihilation picture can be implemented in many well motivated contexts of new physics such as supersymmetry, or extra dimension theories. Phenomenologically, such a scenario can be interesting because the current LHC constraint on the color particle mass is weaker (much less than a TeV) with a compressed spectrum.
The aim of our work is to carry out a precision calculation of the DM relic abundance via coannihilation. To this end, the following three effects must be taken into account.

First, in the annihilation of the colored coannihilating partner and its antiparticle, the gluon exchange between the two initial state colored particles can strongly distort their wavefunction at the origin giving rise to the Sommerfeld effect [2, 3, 4]. During the thermal freeze out, the DM particles and their coannihilating partners in the universe are no longer ultrarelativistic, therefore the Born level cross sections have to be modified to include such a nonperturbative effect. Such an effect has been included previously in [5, 7, 6]. It is worth emphasizing that the Sommerfeld enhancement is in effect above the kinematic threshold of , when there exists a longrange attractive force between two onshell particles.

The second nonperturbative effect lies in the bound state formation channels where the coannihilating partner and its antiparticle radiatively captures each other and form a bound state , before the annihilation takes place. Such a bound state, the QCD onia, exists below the threshold,
(1) which makes it distinct from the Sommerfeld enhancement that is active above , between two onshell particles. Such bound state effects have been found to play an important role in many aspects of DM interacting with a light mediator, from the production in the early universe [8, 9, 10, 11, 12, 14, 13, 15], to indirect detection [16, 17, 18, 19, 20, 21, 22, 23] and direct detection [24, 25, 26], and to collider searches [27, 28, 29, 30]. These motivate us to also include such an effect in the coannihilation calculation.
Based on the Kramers formula [31, 32], generalizable to the QCD case, the most abundantly formed bound states have principle quantum number smaller than , where is the DM relative velocity during the freeze out and numerically it is comparable to . As a result, the calculation including the ground state and first few excited states is expected to capture the dominant contribution.

Third, although the mechanism of coannihilation relies on the DM’s converting into its partners, we also want to emphasize the importance of direct annihilation channels of the DM at low temperature, which are often neglected in the previous simplified model analysis. These channels include the DM selfannihilation and the DMpartner coannihilation channels. Once a DMpartnerquark coupling is introduced for efficient conversion, these channels are automatically induced, and their strength could be fixed in specific models (such as the simple SUSY models discussed below). During the freeze out, the rate of these channels could be competitive to or even more important than the strong interaction channels involving the DM partner. The key reason is that at low enough temperature the DM number density is substantially higher than that of its partner due to the different Boltzmann suppression dictated by their mass difference.
We find that including the annihilation channels involving the DM tends to weaken the relative importance of nonperturbative effects. As a quantitative example, we find that in the binorighthandedsbottom coannihilation case the bound state channel contribution to the overall annihilation rate can be as large as 30%, while in the binorighthandedstop coannihilation case it is always less than 10% because righthanded stop has a larger hypercharge than sbottom.
In addition to the above refinements in the calculation of relic density, we also emphasize the interesting interplay with direct detection of DM and collider searches at the LHC. In particular, we point out that collider searches of stops and sbottoms in the regime of ”compressed spectra,” where the stops and sbottoms are nearly degenerate with the DM particle, invoke the same underlying assumption as the binostop and binosbottom coannihilations. As a result, limits from LHC searches can directly constrain the parameter space in the coannihilation scenarios. This is especially important considering the existence of a ”neutrino floor” in direct detections, which sets a lower limit on the magnitude of direct detection rate that can be probed by conventional experimental setup. For example, we will demonstrate that, in the case of binostop coannihilation, regions of parameter space with a direct detection rate above the neutrino floor is being excluded by the latest limits from the LHC Run 2.
This paper is organized as the following. In section 2, we set up the general formalism for calculating the DM relic abundance in the coannihilation mechanism with all the above effects included. We apply this formalism to calculate the required mass difference between DM and its partner in several simple models including the bino DM coannihilating with a stop or sbottom in section 3, and a vector boson DM coannihilating with a fermionic top partner in section 4. We confront our results with the existing and future experimental probes of the parameter space via DM direct detection and, in the case of binosbottom, the LHC searches. We find including bound state channels could have strong implications in these searches. The conclusion is drawn in section 5.
During the course of this study, several recent works [12, 13, 14] appeared which partially overlap with our discussion on the QCD bound state effects in DM coannihilation. However, in [12, 13], the annihilation channels directly involving the DM (the third bullet in the above) have been neglected in those studies. As will be shown, these neglected effects turn out to be important in the low temperature. On the other hand, [14] considered the class of models where the DMpartner conversion happens through higher dimensional operators, which is a different scenario from the binostop, binosbottom, and vector DM and fermionic top partner coannihilation models considered in this work.
2 The Coannihilation Scenario and Bound State Channels
2.1 The Boltzmann Equations
We first give a general description of the coannihilation scenario and include the bound state effects. We call the DM particle and assume it is a SM singlet and conjugate to itself, and call its heavier partner particle which is colored under the fundamental representation of . In the cases we study, there is a tree level coupling between , and a third generation SM quark, denoted by . Furthermore, the bound states (onia) formed by coannihilation partners and is denoted by . During the DM freeze out, the following processes are relevant

DM selfannihilation into SM (labelled by ):

DMpartner conversion (labelled by ): and , as well as decay and inverse decays of

DMpartner coannihilation (labelled by ):

Partner annihilation with Sommerfeld enhancement (labelled by ): , , etc.

Bound state decay into DM particles (labelled by ):

Bound state decay into SM particles (labelled by ): or , , etc.

Bound state radiative capture and dissociation (labelled by ):
The coupled Boltzmann equations involving are
(2)  
(3)  
(4) 
In the above equations, where is the number density of a species and is the entropy density in the universe. We assume Boltzmann distribution for calculating the thermal number densities. The parameter and is the temperature of the universe. The reaction rate is defined as
(5) 
This expression can be further simplified for the cases of decay and annihilations. See appendix A for more detailed discussion on the form of in decay and annihilation processes. Clearly, is identical for a process and its inverse process.
We assume CP conservation, and , , thus the Boltzmann equation for is the same as that for . This also explains the factor of ’s in Eq. (2), which arises from coannihilating with or converting into the antipartner particle . In contrast, there are no factor of 2 in front of . It is canceled by a symmetry factor for identical particles in the initial state integral, as commented in [1].
The above coupled Boltzmann equations can be further simplified with the following considerations. First, a useful observation is that during freeze out, the bound state formation and decay rates are all much larger than the Hubble expansion rate, i.e., . This implies that it is a good approximation to set the value of equal to the quasistatic solution, which makes the righthand side of Eq. (4) vanishing at every time,
(6) 
Second, after the freeze out, the partner particles will eventually decay into the DM . Therefore, the final DM relic abundance is directly related to the quantity . The other useful observation is that, when the masses of and are close enough, the conversion rate is also much larger than the Hubble expansion rate during freeze out. This implies an approximate relation [1],
(7) 
Using Eqs. (6) and (7), we can add up the Boltzmann equations for and simplify them into the following compact form,
(8) 
where
(9) 
2.2 Bound State Formation and Decays
We first calculate the radiative bound state formation cross section for . Several remarks are in order. First, because the particle is a color triplet, the potential between and is only attractive if they form a color singlet state. This means the bound state must be in a color singlet, and in turn, the initial state must be in a color octet state and they are repulsive to each other. We denote the difference between the initial and final state potential energy as
(10) 
where is the strong coupling for the bound state, evaluated at the energy scale equal to the inverse Bohr radius ( is the principle quantum number) [33], and is the strong coupling for the initial scattering state, evaluated at the energy scale equal to the momentum of the initial particle in the centerofmass frame. In addition, the color factor for radiating a gluon from the or line is 4/3. This could be understood by sandwiching the transition matrix element between the initial state and final state , which gives , where () are fundamental (adjoint) representation color indices. Squaring this yields 4/3. See appendix B for more details. There is an additional factor 1/9 for averaging over the initial state colors. In appendix C, we give a brief discussion of how shows up in the above transition matrix element.
Second, thanks to the nonabelian nature of the gauge group, there is an additional diagram for the capture into bound states, where a real gluon is radiated from the Coulomb gluon exchanged between and . Effectively, it contributes as another potential term
(11) 
Here we evaluate at the energy scale which is the average of the inverse Bohr radius and initial particle momentum. See appendix D for a detailed derivation of this contribution. This piece of contribution was first noticed in [22] and [13].
Under the dipole approximation, the general formation cross section of bound state formation in the force case was derived in [19]. Here, we generalize the result for the case of QCD bound state formation, with arbitrary quantum numbers ,
(12) 
where is the sum of the binding energy of the ’th bound state level and the kinetic energy of incoming state. and are the radial part of the wavefunction of the bound state and initial scattering state, defined as , , respectively. Here is the sum of the potential energies in Eqs. (10) and (11),
(13) 
Because of the contribution from nonabelian gauge interaction, , the above cross section does not agree with the analytic one given in [10, 12] for the case of ground state. In particular, an accidental cancellation occurs when evaluating the above matrix elements for the transition between a low energy scattering state to the ground state, as noticed in [13]. In our calculations below, we take into account of the and bound states. The latter gives about 1020% correction to the former.
This above cross section is derived based on the Schrödinger equation with the spinorbit interactions neglected at leading order, therefore, it is insensitive to the spin of the coannihilating partner and identical between the scalar and fermionic partners. In the latter case, when Eq. (12) is used as the spin averaged cross section for calculating the reaction rate , three quarters of the formed bound states will be in spin triplet () and only one quarter in spin singlet (), i.e., , .
Next, we calculate the bound state decay rates via the annihilation inside it. We define the bound states made of both scalar and fermion colored partners in the appendix E. For wave bound state decay, the general procedure is to first calculate the amplitude of the process final state for at rest with the external legs of amputated ( are their color indices). The bound state decay amplitude to the same final state is obtained by using the projection operators (see also [34]),
where is the bound state wavefucntion at the origin, is the four momentum of the bound state. In the above denote the spin singlet (scalar) state while is the spin triplet (vector) state, and is the polarization vector of the state.
With this approach, we calculate the wave bound state decay rates into gluons, which are (see also [35]),
(14) 
We neglect the electroweak decay channels of the bound states hereafter, which is subdominant and their effects are small. We will present the bound state decay rates into DM particles in specific models in the next sections.
3 BinoStop/Sbottom Coannihilation
In this section, we consider more concrete models where the bino DM (i.e., an and singlet fermion) coannihilates with a righthanded top or bottom scalar quark (which we will call stop or sbottom). This corresponds to the limit in MSSM where all other superpartners are much heavier and decoupled. Notice that this is similar to the assumption of the ”simplified model” approach adopted in most direct searches for stop and sbottom at the LHC. Furthermore, the ”coannihilation” scenario where the bino and stop/sbottom are close in mass corresponds to the ”compressed spectra” scenario [36], where the limits from direct searches are significantly weakened.
3.1 BinoStop
The relevant interacting Lagrangian for the binostop case is
(15) 
where is the hypercharge gauge coupling, is the Higgs doublet and GeV. Applying the general discussion in the previous section to here, we have and . With these interactions, we calculate the freeze out of bino DM in the presence of a slightly heavier stop, by including the bino selfannihilation process , the binostop coannihilation processes , the stopantistop annihilation processes , as well as the stopponium bound state channels. We calculate the Bornlevel cross sections of these processes using CalCHEP [37] which agree with the wave result presented in [6]. We also include the Sommerfeld enhancement/suppression factors discussed in [5, 38].
For the bound state channels, the radiative stopponium formation cross section and its partial decay rate into two gluons have been calculated in the previous section. We consider the ground state in the calculation, whose quantum numbers are , and we call it . There are two additional decay channels, the bound state decay into two binos (via channel top quark exchange) and bound state decay into (via channel bino exchange). Because the bound state is a color singlet, the channel gluon exchange does not contribute to . Their partial decay width are
(16) 
where . For ground state and the strong coupling is evaluated at the energy scale equal to the inverse Bohr radius . Both of the and decay rates contribute to in the Boltzmann equation, and in Eq. (9). We neglected decays into channels because they are always subdominant to .
We solve the Boltzmann equation Eq. (8) with all the above processes included. The parameter space for binostop coannihilation to yield the correct DM relic abundance is depicted in Fig. 1. The black dotdashed curve corresponds to the calculation including Bornlevel stopantistop annihilation and without bound state channels. The blue dotdashed curve corresponds to the calculation including stopantistop annihilation with Sommerfeld effects but without bound state effects. The red solid curves correspond to to the calculation with both Sommerfeld and bound state effects included. Comparing the red and blue curves, we find that including the bound state effects could further increase the stopbino mass difference by 5 GeV for bino mass TeV and as large as more than 10 GeV for TeV.
In our results, we find the over all stopbino mass difference is larger than that found in [12].^{§}^{§}§Our results without the QCD boundstate effects, the dotted and the dashed curves in Fig. 1, agree with those in Ref. [6]. On the other hand, the increase in the mass difference by including bound state effects (compare the red and blue curves) is less significant. A crucial difference we notice is that [12] has only included the stop annihilation channels which are strong interactions but neglected all the electroweak processes, especially those involving the bino. For the binostop coannihilation scenario, we find the electroweak interaction rates to be dominant over the QCD interaction rates during the freeze out for bino mass below 300400 GeV, and this is true for an even wider mass range at late times. This is mainly because the number density of stop (which is heavier) receives much more Boltzmann suppression than that of bino at temperatures smaller than the mass difference. Precision calculation of the coannihilation mechanism requires taking into account of all the strong and electroweak interactions, as well as the nonperturbative Sommerfeld and bound state effects. The importance of the electroweak interaction processes is shown in an example in Fig. 2. Once these processes are included, the bound state effects during freeze out is less than about 10% throughout the bino mass range.
Another bound state effect that has been neglected in the previous analysis is its decay into binos, Eq. (3.1). This channel could have an impact on the behavior of the reaction rate at low temperatures of the universe. When the temperature is large enough (before and around freeze out), we typically find the following hierarchy among the reaction rates in Eq. (9), , and as a result, . In contrast, at lower temperature (around and after freeze out), the hierarchy typically turns into , again because the number density of stop is more suppressed than bino. In this case, . There is a transition in the behavior of between these two regimes. This feature is shown in an example in Fig. 3.
Because of the small mass difference in the stopbino coannihilation region, at the LHC the pairproduced stop has to decay into bino and through an offshell top quark (and offshell boson),
(17) 
The latest LHC constraint on such a decay mode is given in [39], which excludes the bino mass up to GeV for a compressed spectrum. This is shown by the light orange shaded region in Fig. 1. Another potentially useful channel to search for the stop that is nearly degenerate with the bino is through stopponium production at colliders. A recent study shows the future running of high luminosity LHC could reach stop mass up to 350 GeV via the stopponium decay channels [40].
The direct detection processes of bino DM in this case is generated at loop level. The box diagrams could give effective binogluon interactions and the triangle diagrams give binoHiggs interactions. Both will contribute to the isospinconserving and spinindependent scattering on the nucleon target. We take both into account after correcting the sign of the Higgs contribution given in [6]. In Fig. 1, in the region to the left (right) of the green dotted curve, the DMnucleon scattering cross section is above (below) the neutrino floor. It is worth noting that the current LHC stop search has almost ruled out all the parameter space above the neutrino floor. In the future we must resort to the high energy/luminosity colliders for probing the remaining binostop coannihilation parameter space.
3.2 BinoSbottonium
Next, we discuss the case of binosbottom coannihilation. The simplified Lagrangian in this case take the form
(18) 
The binosbottom coannihilation and freeze out is very similar to the binostop case, and early universe calculation results in both cases share many features discussed above, which will not repeated. The parameter space for binosbottom coannihilation to yield the correct DM relic abundance is depicted in Fig. 4. We find the inclusion of the bound state channels has a stronger effect in enlarging the sbottombino mass difference than the stopbino case. As shown in Fig. 5, for low bino mass, the bound state channels could contribute as large as 30% to the overall annihilation rate during freeze out. This is mainly because the righthanded sbottom has smaller hypercharge than the stop, so the electroweak interactions involving the bino have smaller rates.
What we find the most interesting is the phenomenological implication of this result, especially in direct detection. In this case, because the bottom quark is light, the loop diagrams generating binogluon effective interaction (with sbottom and bottom running in the loop) is enhanced when the sbottombino mass difference is small [41, 42, 43, 44]. The enhancement is the strongest when the mass difference is close to the bottom quark mass. This is understood because the virtualness of the sbottom in the loop is controlled by this mass difference. The effective binogluon coupling is only suppressed by the mass square difference instead of the sbottom mass square. In addition, there is no destructive interference like the stop case, because the effective binoHiggs coupling is suppressed by the small bottom Yukawa coupling. In Fig. 4, the green shaded region with low bino mass and small mass difference is ruled out by the current LUX result [45]. The future LZ reach [46] as well as the neutrino floor are shown by the dashed and dotted green curves, respectively. An important message one learns here is that the sbottomium bound state effects pushes part of the coannihilation parameter space (the red solid curve) to lie below the neutrino floor. We will not be able to completely exclude the binosbottom coannihilation scenario by doing the direct detection experiments only.
It is interesting to note that the effects considered in this work give rise to a larger mass difference between the bino and its coannihilating sbottom that would fit the observed relic density, which in turn results in a larger direct detection rate. This is due to the fact that the oneloop diagrams contributing to the gluonbino interactions develops a pole after analytic continuation at [43]
(19) 
The pole exists only when the gluon momentum vanishes, which is a very good approximation in direct detection. Since, by assumption, , the pole is outside of the physical region of the scattering amplitudes. Nevertheless, since , one could be sitting in the vicinity of the pole when the bino and sbottom are almost degenerate in mass and, therefore, receive an enhancement in the direct detection crosssection. Such an enhancement has been studied and demonstrated previously [43, 44]. These considerations makes it clear why the direct detection rate for binosbottom is reduced when the mass difference becomes larger, because one now moves further away from the pole in Eq. (19).
In Fig. 4, we also show the current LHC searches for sbottom decays,
(20) 
which results in the missing energy signature [47], as well as the monojet channel [48, 39] in which case the bottom quarks from sbottom decay are too soft so an initial state jet radiation is triggered on. As emphasized already, the direct searches at the LHC involve the same assumption as the coannihilation scenario, that is the bino is close in mass to its coannihilation partner, while all other superpartners are heavy and decoupled. We see that a bino with mass less than around 600 GeV is excluded by direct searches. With more data collected at the LHC, the direct searches will continue to have an impact on the allowed binosbottom coannihilation region, and might eventually rule out the region where the direct detection rate lies below the neutrino floor.
4 VDMTopPartner Coannihilation
In this section, we consider another scenario with a vector dark matter (VDM) coannihilating with a heavier fermionic top quark partner . This could be motivated by extra dimensional models where the DM stability is related to a KK parity [49]. ^{¶}^{¶}¶Strictly speaking, the spectrum required for the coannihilation studied in this section does not appear in the usual minimal UED models [50]. However, one could potentially introduce branelocalized kinetic terms to generate the desired mass spectrum [51, 52]. We take a simplified model approach and assume the following renormalizable interacting Lagrangian is
(21) 
In the following discussions we will assume to be real and set for simplicity.^{∥}^{∥}∥Modifying this assumption will not alter qualitatively the results shown below. Furthermore, we leave both and as free parameters in the model.
In the coannihilation calculation, we have included for DM self annihilation, for DMpartner coannihilation, and taking into account the Sommerfeld effects from QCD interaction. The channel exchange of particle can also mediate the samesign annihilation of . We have also included the bound state channels. Here because the top partner is a fermion, there are two ground states made of , ^{**}^{**}**It is also possible for to bind into same sign bound state if they form a color sextet [53]. Such a bound state could annihilate decay into two samesign top quarks. Numerically, we find the contribution from such a bound state to coannihilation is negligible, because it has a weaker Coulomb potential () and it has to decay via the weak coupling . spin zero state and spin one state . Their formation cross sections and partial decay rates into gluons have been discussed in Section 2.2. In addition, the bound state could also decay into two vector DM particles (via channel top quark exchange), and into (via channel exchange). These decay rates are
(22) 
Next, for , in the absence of the coupling, the decay is forbidden by parity under the assumption . ^{††}^{††}††The decay channel can be turned on with and both nonzero, which break parity. In practice, we find that it makes very little change to the overall bound state effect, because most of the that are formed are dissociated at much larger rate than decay. The decay rate of is,
(23)  
Our results are shown in Fig. 6, for and , respectively. For , we find that for GeV, the VDM self annihilation cross section alone is large enough to explain the relic abundance, which in turn allows the mass difference between toppartner and VDM to be very large. We find the bound state effect on the toppartnerVDM mass difference is smaller than the squarkbino case. One important reason is that 3/4 of the bound states that are formed in this case are in the state, which unlike , always has to decay into three gluons because of the LandauYang theorem (note is a colorsinglet bound state), which is too slow. As a result, most of the that are formed during the freeze out temperature are quickly dissociated before it has time to decay into gluons. More quantitatively, Eq. (9) tells that the contribution to the bound state rate from is much smaller than .
Next we discuss the interplay of our result with the DM direct detection searches, where the current LUX exclusion, the future LZ reach and the neutrino floor are shown by the green shaded region, green dashed and dotted curves. The spinindependent VDMnucleon scattering could happen via the gluon effective coupling which is generated at loop level [49], and the interaction. ^{‡‡}^{‡‡}‡‡In addition to the tree level term in Eq. (21), there is also a loop level contribution to the coupling which involves the top quark and the top partner. We find the loop contribution is divergent which means that it makes the tree level coupling run with energy scales. In our calculation, the value of we use corresponds to the one after this renormalzation. We vary the value of within the range from 0 to 0.03. The upper limit is motivated by the coupling between Higgs and the first KK hypercharge gauge boson in the minimal UED, which is . We find that for , there is a cancellation between different box diagrams contributing to the DMgluon effective operator (as noted in [54]), which results in a big part of the parameter space (shown by the gray shaded region in Fig. 6) dipping below the neutrino floor of direct detection experiments. Such a cancellation is much less visible when we turn on a large enough value of . In contrast, varying the coupling in the above range has little impact on the coannihilation curves.
On the other hand, at the LHC the only decay channel for the fermionic top partner in our simplified model is
(24) 
where the dark matter particle manifests itself as missing transverse energy (MET). However, as can be seen from Fig. 6, the required mass difference between and to explain the observed relic density is quite small, less than the mass of the boson, which implies can only decay through a highly offshell top into quark and an offshell boson. The resulting signature at the LHC for the pair production ^{x}^{x}xThe top partner needs to be pairproduced because they are odd under KKparity. of is therefore soft jets, soft leptons and a large amount of MET. We are not aware of any public results for searches for fermionic top partners in this particular channel. Needless to say, it would interesting to pursue this channel in future searches.
5 Conclusion
To summarize, in this work, we perform a precision calculation of dark matter relic abundance through the coannihilation mechanism in the context of a few simple models, where the dark matter particle (assumed to be a SM gauge singlet) is accompanied by a slightly heavier and colored top/bottom partner particle. We consider the cases where the quantum numbers of the two dark particles allow a renormalizable interaction with a standard model quark. Given this interaction, the dark matter particles could annihilate during freeze out into standard model particles by: 1) being converted into the partner and then the partners annihilate away through strong interactions, and 2) direct annihilating into quarks through a  and channel partner exchange or coannihilate with a partner particle. In the partnerantipartner annihilation channels belonging to 1), we have include the nonperturbative Sommerfeld as well as the QCD bound state effects in our calculation. In contrast, the annihilation channels belong to 2) are sometimes neglected in previous studies. We point out the importance of these channels — they must be taken into account for a precision calculation of freeze out. Based on these considerations, we derive the cosmologically favored mass difference between the dark matter and its partner.
We confront the derived parameter space for dark matter relic abundance with the present and future experimental searches, which reveals several interesting interplay. For the case of binorighthandedstop coannihilation, most of the parameter space is above the neutrino floor in direct detection has already been excluded by the current LHC stop search. Therefore, we have to rely on the future LHC and high energy colliders instead of direct detection experiments to further probe such a dark matter candidate. On the other hand, for the case of binorighthandedbottom coannihilation, we find the direct detection experiments can serve as powerful probes. However, including the bound state effects in coannihilation pushes part of the parameter space (with dark matter mass between 600 GeV and 1 TeV) to be below the neutrino floor. In this region, the collider searches could play a complementary role.
Our analysis can be implemented to various models with nearby colored partners of the dark matter in a straightforward manner. For example, we include a case study of the vector dark matter with a fermionic top quark partner.
Acknowledgement
This work is supported by the DOE grant DESC0010143 at Northwestern, DOE grant DEFG0212ER41811 at UIC and DOE grant DEAC0206CH11357 at ANL. We acknowledge useful conversations with G. Bodwin and C. Wagner. We also thank the authors of Refs. [13, 14] for communications which led to improvements of our work.
Appendix A Thermal Reaction Rates in Boltzmann Equations
In this appendix, we give the expressions for the thermal averaged reaction rate defined in Eq. (5), for decay and annihilation processes. It is worth noting that is the same for a process and its inverse process.
For a decay process , we have
(25) 
where is the Bessel function of the first kind and is the partial decay rate in the rest frame of .
For an annihilation process ,
(26) 
where and . When the initial state particle masses are equal to each other, i.e., , the above formula can be simplified to
(27) 
where is the relative velocity between and .
Appendix B Color Factors
In this appendix, we discuss several processes of colored partner and interaction, and derive the corresponding color factors in the corresponding amplitudes.
We first define the states. For a system where are the color indices, the color singlet state is obtained by contraction with , while the color octet state is obtained by contraction with . The normalization of these factors are obtained by requiring their square to yield the color multiplicity.
The first process we consider is a photon emission in scattering. Clearly, their color indices do not change in such a process, so the color part of the transition operator is simply . The color factor for color singlet to singlet transition is then . The color factor for color octet to octet transition is then . This is a trivial example.
Second, we consider a gluon emission in scattering. In this case, because the gluon is a color octet, one of the initial and final states must be a color octet and the other is color singlet. The transition operator for gluon radiated from the particle is . And the color factor for octet to singlet transition is . Squaring the amplitude and sum over yields the color factor 4/3, which is used for deriving Eq. (12).
Third, we consider the Coulomb potential between and generated by a gluon exchange. In this case the operator is . For the Coulomb potential between in a color singlet, we have the color factor, . For the Coulomb potential between in a color octet, the color factor is, . This explains why in a color singlet (octet) state are attractive (repulsive) to each other, as well as the normalization of the gluonexchange Coulomb potentials.
Appendix C Capture into Bound States: Gluon Radiated from
In the process ( are the color indices), the effective Hamiltonian for an ultrasoft gluon (with color index ) radiated from or is
(28) 
where is the relative momentum operator between and , and is the gluon field. The matrix element of between the final bound state (color singlet) and the initial scattering state (color octet, carrying adjoint color index ) is then (see the first two diagrams in Fig 7),
(29) 
where , and the color factor in the front is obtained based on the discussions in appendix B. is the polarization vector of the radiated gluon. The wavefunctions and are the eigenstates of the Hamiltonians , and , respectively, with and . Here is the strong coupling for the bound state, evaluated at the energy scale equal to the inverse Bohr radius, and is the strong coupling for the initial scattering state, evaluated at the energy scale equal to the momentum of the initial particle in the centerofmass frame. We define , which yields Eq. (10).
Using , we get
(30) 
Appendix D Capture into Bound States via Nonabelian Gauge Interaction
Because the gluons have their own selfinteractions, it is also possible to radiate the gluon from a gluon propagator exchanged between particles, as shown by the third diagram of Fig 7. In nonrelativistic limit, at leading order, the gluon being exchanged is a soft Coulomb gluon. In this case, the interaction between a radiative gluon and two Coulomb gluons can be written as
(31) 
where is the group structure constant. The corresponding covariant derivative on the fundamental field is, . With these, we calculate the amplitude for the third diagram of Fig 7,
(32) 
where is the momentum running on the gluon propagator.
Going back to the coordinate space, we find this leads to the following interacting Hamiltonian
(33) 
Here the strong coupling arise from the coupling of and with the Coulomb gluons, which is located at the vertex connecting the initial scattering state and the final bound state (next to the color indices “A” and “B” in the third diagram of Fig 7). We decide to evaluate at the energy scale which is the average of the inverse Bohr radius and initial particle momentum.
Sandwiching between the bound and scattering states, we get a new contribution to the transition matrix element
(34) 
Appendix E Bound States and Projection Operators
In this appendix, we define the bound states made of scalar or fermion color particle and its antiparticle, and introduce the corresponding projection operator for calculating the wave bound state decay processes.
For a scalar color particle (stop or sbottom), we define its plane wave state to be . Under this convention, . Then the bound state at rest made of is defined as
(36) 
where are the color indices, and is the Fourier transformation of the bound state wavefunction in the momentum space, . It is straightforward to verify that the inner product of two bound state at rest is,