A Reappraisal on Dark Matter Co-annihilating with a Top/Bottom Partner

TL;DR

This paper revisits dark matter relic density calculations involving co-annihilation with top or bottom partners, emphasizing the importance of QCD bound-state effects and electroweak processes, and compares predictions with experimental constraints.

Contribution

It provides a comprehensive framework including bound-state and electroweak effects for co-annihilation scenarios, improving relic density predictions and confronting them with current experimental data.

Findings

Bound-state effects significantly impact relic density calculations.

LHC and direct detection constraints exclude some parameter regions.

Electroweak processes are crucial in low mass co-annihilation scenarios.

Abstract

We revisit the calculation of relic density of dark matter particles co-annihilating with a top or bottom partner, by properly including the QCD bound-states (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: bino-stop and bino-sbottom co-annihilations in supersymmetry, and a vector dark matter co-annihilating 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 bino-stop and bino-sbottom scenarios. In particular, in the bino-stop case recent LHC limits have excluded regions of parameter with a direct detection rate that is above the neutrino floor.

Figures

Equations

$$m_B<2m_Y,$$

$$sHz\frac{d{\mathcal{Y}}_X}{dz}=-{\gamma }_{XX}\left[{\left(\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\right)}^2-1\right]-2{\gamma }_{XY}\left[\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}-1\right]-2{\gamma }_{X\leftrightarrow Y}\left[\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}-\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}\right]-2{\gamma }_{B\leftrightarrow X}\left[{\left(\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\right)}^2-\frac{{\mathcal{Y}}_B}{{\mathcal{Y}}_B^{eq}}\right],$$

$$sHz\frac{d{\mathcal{Y}}_Y}{dz}=-{\gamma }_{Y\bar{Y}}\left[{\left(\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}\right)}^2-1\right]-{\gamma }_{XY}\left[\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}-1\right]+{\gamma }_{X\leftrightarrow Y}\left[\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}-\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}\right]-{\gamma }_{B\leftrightarrow Y}\left[{\left(\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}\right)}^2-\frac{{\mathcal{Y}}_B}{{\mathcal{Y}}_B^{eq}}\right],$$

$$sHz\frac{d{\mathcal{Y}}_B}{dz}=-{\gamma }_{B\leftrightarrow X}\left[\frac{{\mathcal{Y}}_B}{{\mathcal{Y}}_B^{eq}}-{\left(\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\right)}^2\right]-{\gamma }_{B\leftrightarrow Y}\left[\frac{{\mathcal{Y}}_B}{{\mathcal{Y}}_B^{eq}}-{\left(\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}\right)}^2\right]-{\gamma }_{B\leftrightarrow \mathrm{SM}}\left[\frac{{\mathcal{Y}}_B}{{\mathcal{Y}}_B^{eq}}-1\right].$$

$$\gamma$$

$$\frac{{\mathcal{Y}}_B}{{\mathcal{Y}}_B^{eq}}\simeq \frac{{\gamma }_{B\leftrightarrow X}{\left(\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\right)}^2+{\gamma }_{B\leftrightarrow Y}{\left(\frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}\right)}^2+{\gamma }_{B\leftrightarrow \mathrm{SM}}}{{\gamma }_{B\leftrightarrow Y}+{\gamma }_{B\leftrightarrow X}+{\gamma }_{B\leftrightarrow \mathrm{SM}}}.$$

$$\frac{{\mathcal{Y}}_X}{{\mathcal{Y}}_X^{eq}}\simeq \frac{{\mathcal{Y}}_Y}{{\mathcal{Y}}_Y^{eq}}.$$

$$sHz\frac{d{\mathcal{Y}}_{\mathrm{dark}}}{dz}=-\left[{\gamma }_{XX}+4{\gamma }_{XY}+2{\gamma }_{Y\bar{Y}}+2{\gamma }_{\mathrm{bound}\mathrm{state}}\right]\left[{\left(\frac{{\mathcal{Y}}_{\mathrm{dark}}}{{\mathcal{Y}}_{\mathrm{dark}}^{eq}}\right)}^2-1\right],$$

$${\gamma }_{\mathrm{bound}\mathrm{state}}=\sum _B\frac{({\gamma }_{B\leftrightarrow X}+{\gamma }_{B\leftrightarrow Y}){\gamma }_{B\leftrightarrow \mathrm{SM}}}{{\gamma }_{B\leftrightarrow X}+{\gamma }_{B\leftrightarrow Y}+{\gamma }_{B\leftrightarrow \mathrm{SM}}}.$$

$$\Delta V_1(r)=\frac{4{\alpha }_S^{(f)}}{3r}+\frac{{\alpha }_S^{(i)}}{6r},$$

$$\Delta V_2(r)=\frac{3{\alpha }_S}{2r}.$$

$$(\sigma v_{\mathrm{rel}}{)}_{Y\bar{Y}\to Bg}^{n,\ell }=\frac{4{\alpha }_D}{81\pi }{\omega }_n\left[\ell {\left|\int dr{r}^3\left({\omega }_n-\Delta V\right)R_{n\ell }{R}_{k\ell -1}\right|}^2+(\ell +1){\left|\int dr{r}^3\left({\omega }_n-\Delta V\right)R_{n\ell }{R}_{k\ell +1}\right|}^2\right],$$

$$\Delta V(r)=\Delta V_1(r)+\Delta V_2(r).$$

$${\mathcal{A}}_{B\to \mathrm{final}}=\{\begin{array}{ll}\sqrt{\frac{1}{3m_Y}}\Psi (0){\delta }_{ij}{\mathcal{A}}_{ij},& Y\mathrm{is}\mathrm{scalar}\\ \sqrt{\frac{1}{24m_Y^3}}\Psi (0){\delta }_{ij}\mathrm{Tr}\left[\left(\frac{/P}{2}+m_Y\right){\gamma }_5\left(\frac{/P}{2}-m_Y\right){\mathcal{A}}_{ij}\right],& Y\mathrm{is}\mathrm{fermion},B={\eta }_Y\\ \sqrt{\frac{1}{24m_Y^3}}\Psi (0){\delta }_{ij}\mathrm{Tr}\left[\left(\frac{/P}{2}+m_Y\right)/\epsilon (P)\left(\frac{/P}{2}-m_Y\right){\mathcal{A}}_{ij}\right],& Y\mathrm{is}\mathrm{fermion},B={{\rm Y}}_Y\end{array}$$

$${\Gamma }_{B\to gg}=\{\begin{array}{ll}\frac{4\pi {\alpha }_S(m_Y{)}^2}{3m_Y^2}|\Psi (0){|}^2,& Y\mathrm{is}\mathrm{scalar}\\ \frac{8\pi {\alpha }_S(m_Y{)}^2}{3m_Y^2}|\Psi (0){|}^2,& Y\mathrm{is}\mathrm{fermion},B={\eta }_Y\end{array}$$

$${\Gamma }_{B\to ggg}=\frac{40({\pi }^2-9){\alpha }_S(m_Y{)}^3}{81m_Y^2}|\Psi (0){|}^2,Y\mathrm{is}\mathrm{fermion},B={{\rm Y}}_Y$$

$${\mathcal{L}}_{\mathrm{bino}-\mathrm{stop}}=\left(\frac{\sqrt{2}{g}_1}{3}\bar{t}(1-{\gamma }_5)\chi \tilde{t}+h.c.\right)+\frac{2m_t^2}{v^2}{H}^{†}H{\tilde{t}}^{\ast }\tilde{t},$$

$${\Gamma }_{B_{\tilde{t}}\to \chi \chi }=\frac{128\pi {\alpha }_1^2{m}_{\chi }^2{\left(m_{\tilde{t}}^2-m_{\chi }^2\right)}^{3/2}}{27m_{\tilde{t}}^3{\left(m_{\tilde{t}}^2-m_{\chi }^2+m_t^2\right)}^2}|\Psi (0){|}^2,$$

$${\Gamma }_{B_{\tilde{t}}\to t\bar{t}}=\frac{64\pi {\alpha }_1^2{m}_t^2{\left(m_{\tilde{t}}^2-m_t^2\right)}^{3/2}}{81m_{\tilde{t}}^3{\left(m_{\tilde{t}}^2-m_t^2+m_{\chi }^2\right)}^2}|\Psi (0){|}^2.$$

$$\tilde{t}\to t^{\ast }+\chi \to W^{\ast }+b+\chi ,$$

$${\mathcal{L}}_{\mathrm{bino}-\mathrm{sbottom}}=\left(-\frac{\sqrt{2}{g}_1}{6}\bar{b}(1-{\gamma }_5)\chi \tilde{b}+h.c.\right)+\frac{2m_b^2}{v^2}{H}^{†}H{\tilde{b}}^{\ast }\tilde{b}.$$

$$m_{\tilde{b}}+m_b=m_{\chi }.$$

$$\tilde{b}\to b+\chi ,$$

$${\mathcal{L}}_{\mathrm{VDM}-\mathrm{TP}}=\left(\bar{T}{\gamma }^{\mu }\left({\lambda }_{VTt}+{\lambda }_{VTt}^{\prime }{\gamma }_5\right)t{V}_{\mu }+h.c.\right)+{\lambda }_{VVHH}{V}_{\mu }{V}^{\mu }{H}^{†}H.$$

$${\Gamma }_{{\eta }_T\to VV}$$

$${\Gamma }_{{\eta }_T\to t\bar{t}}$$

$${\Gamma }_{{{\rm Y}}_T\to t\bar{t}}$$

$$T\to t+V,$$

$${\gamma }_{a\leftrightarrow cd}=\frac{T{m}_a^2}{2{\pi }^2}{K}_1\left(\frac{m_a}{T}\right)g_a{\Gamma }_{a\to cd},$$

$${\gamma }_{ab\leftrightarrow cd}=\frac{T}{8{\pi }^4}{\int }_{s_{\mathrm{th}}}^{\infty }ds\sqrt{s}{K}_1\left(\frac{\sqrt{s}}{T}\right)|{\mathbf{p}}_a^{\mathrm{cm}}{|}^2{g}_a{g}_b{\sigma }_{a+b\to c+d},$$

$${\gamma }_{ab\leftrightarrow cd}=\frac{T}{64{\pi }^4}{\int }_{s_{\mathrm{th}}}^{\infty }dss\sqrt{s-4m_a^2}{K}_1\left(\frac{\sqrt{s}}{T}\right)\times g_a{g}_b(\sigma v_{\mathrm{rel}}{)}_{ab\to cd},$$

$$H_{int}^{(1)}=\frac{g_S{T}_{i^{\prime }i}^C{\delta }_{j^{\prime }j}}{m_Y}\vec{k}\cdot {\vec{G}}^C(\vec{r}/2)+\frac{g_S{\delta }_{i^{\prime }i}{T}_{j^{\prime }j}^C}{m_Y}\vec{k}\cdot {\vec{G}}^C(-\vec{r}/2)$$

$$V_{fi}^{(1)}=\frac{{\delta }^{CD}}{\sqrt{6}}\frac{g_S}{\mu }⟨{\Psi }_f|\vec{k}|{\Psi }_i⟩\cdot {\vec{\epsilon }}_{\lambda }^C,$$

Full text

A Reappraisal on Dark Matter Co-annihilating

with a Top/Bottom Partner

**Wai-Yee Keung 1 *keung@uic.edu,   Ian Low 2,3 †††ilow@northwestern.edu,   Yue Zhang 2 ‡‡‡yuezhang@northwestern.edu

1Physics Department, University of Illinois at Chicago, IL 60607, USA

2Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA

3High Energy Physics Division, Argonne National Laboratory, Argonne, IL 60439, USA

Abstract: We revisit the calculation of relic density of dark matter particles co-annihilating with a top or bottom partner, by properly including the QCD bound-states (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: bino-stop and bino-sbottom co-annihilations in supersymmetry, and a vector dark matter co-annihilating 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 bino-stop and bino-sbottom scenarios. In particular, in the bino-stop case recent LHC limits have excluded regions of parameter with a direct detection rate that is above the neutrino floor.

NUHEP-TH/17-02

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 $X$ is accompanied in the mass spectrum with a heavier partner $Y$ which is colored and participates in the strong interaction. Both are odd under a $Z_{2}$ 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 so-called co-annihilation 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 DM-partner mass difference (typically less than $\sim$ 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 co-annihilation 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 co-annihilation. To this end, the following three effects must be taken into account.

• •

First, in the annihilation of the colored co-annihilating partner $Y$ 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 co-annihilating partners in the universe are no longer ultra-relativistic, therefore the Born level cross sections have to be modified to include such a non-perturbative 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 $2m_{Y}$, when there exists a long-range attractive force between two on-shell $Y$ particles.

• •

The second non-perturbative effect lies in the bound state formation channels where the co-annihilating partner $Y$ and its anti-particle radiatively captures each other and form a bound state $B$, before the annihilation takes place. Such a bound state, the QCD onia, exists below the $2m_{Y}$ threshold,

[TABLE]

which makes it distinct from the Sommerfeld enhancement that is active above $2m_{Y}$, between two on-shell $Y$ 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 co-annihilation 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 $\sim$ $\alpha_{S}/v$, where $v$ is the DM relative velocity during the freeze out and numerically it is comparable to $\alpha_{S}$. 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 co-annihilation 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 self-annihilation and the DM-partner co-annihilation channels. Once a DM-partner-quark 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 non-perturbative effects. As a quantitative example, we find that in the bino-right-handed-sbottom co-annihilation case the bound state channel contribution to the overall annihilation rate can be as large as 30%, while in the bino-right-handed-stop co-annihilation case it is always less than 10% because right-handed 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 bino-stop and bino-sbottom co-annihilations. As a result, limits from LHC searches can directly constrain the parameter space in the co-annihilation 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 bino-stop co-annihilation, 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 co-annihilation 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 co-annihilating with a stop or sbottom in section 3, and a vector boson DM co-annihilating 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 bino-sbottom, 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 co-annihilation. 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 DM-partner conversion happens through higher dimensional operators, which is a different scenario from the bino-stop, bino-sbottom, and vector DM and fermionic top partner co-annihilation models considered in this work.

2 The Co-annihilation Scenario and Bound State Channels

2.1 The Boltzmann Equations

We first give a general description of the co-annihilation scenario and include the bound state effects. We call the DM particle $X$ and assume it is a SM singlet and conjugate to itself, and call its heavier partner particle $Y$ which is colored under the fundamental representation of $SU(3)_{c}$. In the cases we study, there is a tree level coupling between $X$, $Y$ and a third generation SM quark, denoted by $q$. Furthermore, the bound states (onia) formed by co-annihilation partners $Y$ and $\bar{Y}$ is denoted by $B$. During the DM freeze out, the following processes are relevant

• •

DM self-annihilation into SM (labelled by $XX$): $XX\leftrightarrow q\bar{q}$

• •

DM-partner conversion (labelled by $X\leftrightarrow Y$): $Xq\leftrightarrow Yg$ and $Xg\leftrightarrow Y\bar{q}$, as well as decay and inverse decays of $Y\to X+{\rm SM\ particles}$

• •

DM-partner co-annihilation (labelled by $XY$): $XY\leftrightarrow qg$

• •

Partner annihilation with Sommerfeld enhancement (labelled by $Y\bar{Y}$): $Y\bar{Y}\leftrightarrow gg$, $Y\bar{Y}\leftrightarrow q\bar{q}$, etc.

• •

Bound state $B$ decay into DM particles (labelled by $B\leftrightarrow X$): $B\leftrightarrow XX$

• •

Bound state $B$ decay into SM particles (labelled by $B\leftrightarrow{\rm SM}$): $B\leftrightarrow 2g$ or $3g$, $B\leftrightarrow q\bar{q}$, etc.

• •

Bound state $B$ radiative capture and dissociation (labelled by $B\leftrightarrow Y$): $Y\bar{Y}\leftrightarrow Bg$

The coupled Boltzmann equations involving $X,Y,B$ are

[TABLE]

In the above equations, $\mathcal{Y}\equiv n/s$ where $n$ is the number density of a species and $s$ is the entropy density in the universe. We assume Boltzmann distribution for calculating the thermal number densities. The parameter $z\equiv m_{X}/T$ and $T$ is the temperature of the universe. The reaction rate $\gamma$ is defined as

[TABLE]

This expression can be further simplified for the cases of decay and $2\to 2$ annihilations. See appendix A for more detailed discussion on the form of $\gamma$ in decay and annihilation processes. Clearly, $\gamma$ is identical for a process and its inverse process.

We assume CP conservation, $\mathcal{Y}_{Y}=\mathcal{Y}_{\bar{Y}}$ and $\gamma_{X\bar{Y}}=\gamma_{XY}$, $\gamma_{X\leftrightarrow\bar{Y}}=\gamma_{X\leftrightarrow Y}$, thus the Boltzmann equation for $\mathcal{Y}_{\bar{Y}}$ is the same as that for $\mathcal{Y}_{Y}$. This also explains the factor of $2$’s in Eq. (2), which arises from $X$ co-annihilating with or converting into the anti-partner particle $\bar{Y}$. In contrast, there are no factor of 2 in front of $\gamma_{XX}$. It is canceled by a symmetry factor $1/2$ for identical $X$ 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., $\gamma_{B\leftrightarrow X,Y,g}\gg sHz$. This implies that it is a good approximation to set the value of $\mathcal{Y}_{B}$ equal to the quasi-static solution, which makes the right-hand side of Eq. (4) vanishing at every time,

[TABLE]

Second, after the freeze out, the partner particles $Y,\bar{Y}$ will eventually decay into the DM $X$. Therefore, the final DM relic abundance is directly related to the quantity $\mathcal{Y}_{\rm dark}=\mathcal{Y}_{X}+\mathcal{Y}_{Y}+\mathcal{Y}_{\bar{Y}}$. The other useful observation is that, when the masses of $X$ and $Y$ are close enough, the $X\leftrightarrow Y$ conversion rate is also much larger than the Hubble expansion rate during freeze out. This implies an approximate relation [1],

[TABLE]

Using Eqs. (6) and (7), we can add up the Boltzmann equations for $\mathcal{Y}_{X},\mathcal{Y}_{Y},\mathcal{Y}_{\bar{Y}}$ and simplify them into the following compact form,

[TABLE]

where

[TABLE]

2.2 Bound State Formation and Decays

We first calculate the radiative bound state formation cross section for $Y\bar{Y}\to Bg$. Several remarks are in order. First, because the $Y$ particle is a color triplet, the potential between $Y$ and $\bar{Y}$ is only attractive if they form a color singlet state. This means the bound state $B$ must be in a color singlet, and in turn, the initial state $Y\bar{Y}$ 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

[TABLE]

where $\alpha_{S}^{(f)}$ is the strong coupling for the bound state, evaluated at the energy scale equal to the inverse Bohr radius $\mu=a_{n}^{-1}=2\alpha_{S}m_{Y}/(3n)$ ($n$ is the principle quantum number) [33], and $\alpha_{S}^{(i)}$ is the strong coupling for the initial scattering state, evaluated at the energy scale equal to the momentum of the initial $Y$ particle in the center-of-mass frame. In addition, the color factor for radiating a gluon from the $Y$ or $\bar{Y}$ line is 4/3. This could be understood by sandwiching the transition matrix element $(T^{A})_{ik}\delta_{jl}$ between the initial state $\sqrt{2}(T^{B})_{kl}$ and final state $\delta_{ij}/\sqrt{3}$, which gives $\delta^{AB}/\sqrt{6}$, where $i,j,k,l$ ($A,B$) 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 $\Delta V_{1}$ shows up in the above transition matrix element.

Second, thanks to the non-abelian nature of the $SU(3)_{c}$ 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 $Y$ and $\bar{Y}$. Effectively, it contributes as another potential term

[TABLE]

Here we evaluate $\alpha_{s}$ at the energy scale which is the average of the inverse Bohr radius and initial $Y$ 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 $U(1)$ force case was derived in [19]. Here, we generalize the result for the case of QCD bound state formation, with arbitrary quantum numbers $n,\ell$,

[TABLE]

where $\omega_{n}=E_{n}+k^{2}/(2\mu)$ is the sum of the binding energy of the $n$’th bound state level and the kinetic energy of incoming state. $R_{n\ell}$ and $R_{k\ell}$ are the radial part of the wavefunction of the bound state and initial scattering state, defined as $\Psi_{n}({\bf r})=\sum_{\ell m}R_{n\ell}(r)Y_{\ell m}(\hat{r})$, $\Psi_{k}({\bf r})=\sum_{\ell m}R_{k\ell}(r)Y_{\ell m}(\hat{r})Y^{*}_{\ell m}(\hat{k})$, respectively. Here $\Delta V$ is the sum of the potential energies in Eqs. (10) and (11),

[TABLE]

Because of the contribution from non-abelian gauge interaction, $\Delta V_{2}$, 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 $1S$ and $2S$ bound states. The latter gives about 10-20% correction to the former.

This above cross section is derived based on the Schrödinger equation with the spin-orbit interactions neglected at leading order, therefore, it is insensitive to the spin of the co-annihilating 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 $\gamma_{Y\bar{Y}\to Bg}$, three quarters of the formed bound states will be in spin triplet ($\Upsilon_{Y}$) and only one quarter in spin singlet ($\eta_{Y}$), i.e., $\gamma_{Y\bar{Y}\to\Upsilon_{Y}g}=\frac{3}{4}\gamma_{Y\bar{Y}\to Bg}$, $\gamma_{Y\bar{Y}\to\eta_{Y}g}=\frac{1}{4}\gamma_{Y\bar{Y}\to Bg}$.

Next, we calculate the bound state decay rates via the $Y\bar{Y}$ annihilation inside it. We define the bound states made of both scalar and fermion colored partners in the appendix E. For $S$-wave bound state decay, the general procedure is to first calculate the amplitude $\mathcal{A}_{ij}$ of the process $Y\bar{Y}\to$ final state for $Y,\bar{Y}$ at rest with the external legs of $Y,\bar{Y}$ amputated ($i,j$ are their color indices). The bound state decay amplitude to the same final state is obtained by using the projection operators (see also [34]),

[TABLE]

where $\Psi(0)$ is the bound state wavefucntion at the origin, $P^{\mu}$ is the four momentum of the bound state. In the above $\eta_{Y}$ denote the spin singlet (scalar) state while $\Upsilon_{Y}$ is the spin triplet (vector) state, and $\varepsilon$ is the polarization vector of the $\Upsilon_{Y}$ state.

With this approach, we calculate the $S$-wave bound state decay rates into gluons, which are (see also [35]),

[TABLE]

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 $X$ in specific models in the next sections.

3 Bino-Stop/Sbottom Co-annihilation

In this section, we consider more concrete models where the bino DM (i.e., an $SU(3)_{c}$ and $SU(2)_{L}$ singlet fermion) $\chi$ co-annihilates with a right-handed 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 ”co-annihilation” 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 Bino-Stop

The relevant interacting Lagrangian for the bino-stop case is

[TABLE]

where $g_{1}$ is the hypercharge gauge coupling, $H$ is the Higgs doublet and $v=246\,$GeV. Applying the general discussion in the previous section to here, we have $X=\chi$ and $Y=\tilde{t}$. With these interactions, we calculate the freeze out of bino DM in the presence of a slightly heavier stop, by including the bino self-annihilation process $\chi\chi\to t\bar{t}$, the bino-stop co-annihilation processes $\chi\tilde{t}\to tg,th,tZ,t\gamma,bW$, the stop-anti-stop annihilation processes $\tilde{t}^{*}\tilde{t}\to gg,hh,ZZ,WW,\gamma\gamma,gh,gZ,g\gamma,hZ,g\gamma,Z\gamma$, as well as the stopponium bound state channels. We calculate the Born-level cross sections of these processes using CalCHEP [37] which agree with the $S$-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 $J^{PC}=0^{++}$, and we call it $B=B_{\tilde{t}}$. There are two additional decay channels, the bound state decay into two binos (via $t$-channel top quark exchange) and bound state decay into $t\bar{t}$ (via $t$-channel bino exchange). Because the bound state is a color singlet, the $s$-channel gluon exchange does not contribute to $B\to t\bar{t}$. Their partial decay width are

[TABLE]

where $\alpha_{1}=g_{1}^{2}/(4\pi)$. For ground state $|\Psi(0)|^{2}=8\alpha_{S}m_{\tilde{t}}^{3}/(27\pi)$ and the strong coupling $\alpha_{S}$ is evaluated at the energy scale equal to the inverse Bohr radius $\mu=2\alpha_{S}m_{Y}/3$. Both of the $B_{\tilde{t}}\to t\bar{t}$ and $B_{\tilde{t}}\to gg$ decay rates contribute to $\gamma_{B\leftrightarrow{\rm SM}}$ in the Boltzmann equation, and in Eq. (9). We neglected $B_{\tilde{t}}$ decays into $Z,\gamma,h$ channels because they are always subdominant to $B_{\tilde{t}}\to gg$.

We solve the Boltzmann equation Eq. (8) with all the above processes included. The parameter space for bino-stop co-annihilation to yield the correct DM relic abundance is depicted in Fig. 1. The black dot-dashed curve corresponds to the calculation including Born-level stop-anti-stop annihilation and without bound state channels. The blue dot-dashed curve corresponds to the calculation including stop-anti-stop 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 stop-bino mass difference by 5 GeV for bino mass $m_{\chi}\lesssim 1\,$TeV and as large as more than 10 GeV for $m_{\chi}\sim 2\,$TeV.

In our results, we find the over all stop-bino mass difference is larger than that found in [12].§§§Our results without the QCD bound-state 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 bino-stop co-annihilation scenario, we find the electroweak interaction rates to be dominant over the QCD interaction rates during the freeze out for bino mass below $\sim$ 300-400 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 co-annihilation mechanism requires taking into account of all the strong and electroweak interactions, as well as the non-perturbative 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 $\gamma_{\rm bound\ state}$ 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), $\gamma_{B\leftrightarrow\tilde{t}}\gg\gamma_{B\leftrightarrow{\rm SM}}\gg\gamma_{B\leftrightarrow\chi}$, and as a result, $\gamma_{\rm bound\ state}\simeq\gamma_{B\leftrightarrow{\rm SM}}$. In contrast, at lower temperature (around and after freeze out), the hierarchy typically turns into $\gamma_{B\leftrightarrow{\rm SM}}\gg\gamma_{B\leftrightarrow\chi}\gg\gamma_{B\leftrightarrow\tilde{t}}$, again because the number density of stop is more suppressed than bino. In this case, $\gamma_{\rm bound\ state}\simeq\gamma_{B\leftrightarrow\chi}$. There is a transition in the behavior of $\gamma_{\rm bound\ state}$ between these two regimes. This feature is shown in an example in Fig. 3.

Because of the small mass difference in the stop-bino co-annihilation region, at the LHC the pair-produced stop has to decay into bino and through an off-shell top quark (and off-shell $W$-boson),

[TABLE]

The latest LHC constraint on such a decay mode is given in [39], which excludes the bino mass up to $\sim 500\,$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 $\to\gamma\gamma,Z\gamma$ decay channels [40].

The direct detection processes of bino DM in this case is generated at loop level. The box diagrams could give effective bino-gluon interactions and the triangle diagrams give bino-Higgs interactions. Both will contribute to the isospin-conserving and spin-independent 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 DM-nucleon 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 bino-stop co-annihilation parameter space.

3.2 Bino-Sbottonium

Next, we discuss the case of bino-sbottom co-annihilation. The simplified Lagrangian in this case take the form

[TABLE]

The bino-sbottom co-annihilation and freeze out is very similar to the bino-stop case, and early universe calculation results in both cases share many features discussed above, which will not repeated. The parameter space for bino-sbottom co-annihilation 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 sbottom-bino mass difference than the stop-bino 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 right-handed 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 bino-gluon effective interaction (with sbottom and bottom running in the loop) is enhanced when the sbottom-bino 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 bino-gluon 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 bino-Higgs 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 co-annihilation parameter space (the red solid curve) to lie below the neutrino floor. We will not be able to completely exclude the bino-sbottom co-annihilation 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 co-annihilating 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 one-loop diagrams contributing to the gluon-bino interactions develops a pole after analytic continuation at [43]

[TABLE]

The pole exists only when the gluon momentum vanishes, which is a very good approximation in direct detection. Since, by assumption, $m_{\chi}<m_{\tilde{b}}$, the pole is outside of the physical region of the scattering amplitudes. Nevertheless, since $m_{b}\ll m_{\chi},m_{\tilde{b}}$, 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 cross-section. Such an enhancement has been studied and demonstrated previously [43, 44]. These considerations makes it clear why the direct detection rate for bino-sbottom 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,

[TABLE]

which results in the $2b+$ 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 co-annihilation scenario, that is the bino is close in mass to its co-annihilation 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 bino-sbottom co-annihilation region, and might eventually rule out the region where the direct detection rate lies below the neutrino floor.

4 VDM-Top-Partner Co-annihilation

In this section, we consider another scenario with a vector dark matter (VDM) $V_{\mu}$ co-annihilating with a heavier fermionic top quark partner $T$. 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 co-annihilation studied in this section does not appear in the usual minimal UED models [50]. However, one could potentially introduce brane-localized kinetic terms to generate the desired mass spectrum [51, 52]. We take a simplified model approach and assume the following renormalizable interacting Lagrangian is

[TABLE]

In the following discussions we will assume $\lambda_{VTt}$ to be real and set $\lambda^{\prime}_{VTt}=0$ for simplicity.∥∥∥Modifying this assumption will not alter qualitatively the results shown below. Furthermore, we leave both $\lambda_{VTt}$ and $\lambda_{VVHH}$ as free parameters in the model.

In the co-annihilation calculation, we have included $VV\to t\bar{t},hh$ for DM self annihilation, $VT\to gt$ for DM-partner coannihilation, and $T\bar{T}\to gg,q\bar{q}$ taking into account the Sommerfeld effects from QCD interaction. The $t$-channel exchange of $V$ particle can also mediate the same-sign annihilation of $TT\to tt$. We have also included the bound state channels. Here because the top partner $T$ is a fermion, there are two ground states made of $T\bar{T}$, ******It is also possible for $TT$ to bind into same sign bound state if they form a color sextet [53]. Such a bound state could annihilate decay into two same-sign top quarks. Numerically, we find the contribution from such a bound state to co-annihilation is negligible, because it has a weaker Coulomb potential ($V=-\alpha_{S}/(3r)$) and it has to decay via the weak coupling $\lambda_{VTt}$. spin zero state $\eta_{T}$ and spin one state $\Upsilon_{T}$. 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 $t$-channel top quark exchange), and into $t\bar{t}$ (via $t$-channel $V$ exchange). These decay rates are

[TABLE]

Next, for $\Upsilon_{T}$, in the absence of the $\lambda^{\prime}_{VTt}$ coupling, the decay $\Upsilon_{T}\to VV$ is forbidden by $C$-parity under the assumption $\lambda^{\prime}_{VTt}=0$. ††††††The $\Upsilon_{T}\to VV$ decay channel can be turned on with $\lambda_{VTt}$ and $\lambda^{\prime}_{VTt}$ both nonzero, which break $C$-parity. In practice, we find that it makes very little change to the overall bound state effect, because most of the $\Upsilon_{T}$ that are formed are dissociated at much larger rate than decay. The decay rate of $\Upsilon_{T}\to t\bar{t}$ is,

[TABLE]

Our results are shown in Fig. 6, for $\lambda_{VTt}=\{0.3,0.6\}$ and $\lambda_{VVHH}=\{0,0.01,0.03\}$, respectively. For $\lambda_{VTt}=-0.6$, we find that for $m_{V}\lesssim 400\,$GeV, the VDM self annihilation cross section alone is large enough to explain the relic abundance, which in turn allows the mass difference between top-partner and VDM to be very large. We find the bound state effect on the top-partner-VDM mass difference is smaller than the squark-bino case. One important reason is that 3/4 of the bound states that are formed in this case are in the $\Upsilon_{T}$ state, which unlike $\eta_{T}$, always has to decay into three gluons because of the Landau-Yang theorem (note $\Upsilon_{T}$ is a color-singlet bound state), which is too slow. As a result, most of the $\Upsilon_{T}$ 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 $\Upsilon_{T}$ is much smaller than $\eta_{T}$.

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 spin-independent VDM-nucleon scattering could happen via the $V$-gluon effective coupling which is generated at loop level [49], and the $VVHH$ interaction. ‡‡‡‡‡‡In addition to the tree level $\lambda_{VVHH}$ term in Eq. (21), there is also a loop level contribution to the $VVHH$ 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 $\lambda_{VVHH}$ we use corresponds to the one after this renormalzation. We vary the value of $\lambda_{VVHH}$ 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 $g_{1}^{2}/4\simeq 0.03$. We find that for $\lambda_{VVHH}=0$, there is a cancellation between different box diagrams contributing to the DM-gluon 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 $\lambda_{VVHH}$. In contrast, varying the $\lambda_{VVHH}$ coupling in the above range has little impact on the co-annihilation curves.

On the other hand, at the LHC the only decay channel for the fermionic top partner in our simplified model is

[TABLE]

where the dark matter particle $V$ manifests itself as missing transverse energy (MET). However, as can be seen from Fig. 6, the required mass difference between $V$ and $T$ to explain the observed relic density is quite small, less than the mass of the $W$ boson, which implies $T$ can only decay through a highly off-shell top into $b$-quark and an off-shell $W$ boson. The resulting signature at the LHC for the pair production xxxThe top partner needs to be pair-produced because they are odd under KK-parity. of $T$ 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 co-annihilation 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 $t$- and $u$-channel partner exchange or co-annihilate with a partner particle. In the partner-anti-partner annihilation channels belonging to 1), we have include the non-perturbative 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 bino-right-handed-stop co-annihilation, 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 bino-right-handed-bottom co-annihilation, we find the direct detection experiments can serve as powerful probes. However, including the bound state effects in co-annihilation 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 DE-SC0010143 at Northwestern, DOE grant DE-FG-02-12ER41811 at UIC and DOE grant DE-AC02-06CH11357 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 $\gamma$ defined in Eq. (5), for $1\to 2$ decay and $2\to 2$ annihilation processes. It is worth noting that $\gamma$ is the same for a process and its inverse process.

For a decay process $a\to c+d$, we have

[TABLE]

where $K_{1}$ is the Bessel function of the first kind and $\Gamma_{a\to cd}$ is the partial decay rate in the rest frame of $a$.

For an annihilation process $a+b\to c+d$,

[TABLE]

where $s_{\rm th}={\rm Max}[(m_{a}+m_{b})^{2},(m_{c}+m_{d})^{2}]$ and $|{\bf p}_{a}^{\rm cm}|=\sqrt{[s-(m_{a}+m_{b})^{2}][s-(m_{a}-m_{b})^{2}]}/(2\sqrt{s})$. When the initial state particle masses are equal to each other, i.e., $m_{a}=m_{b}$, the above formula can be simplified to

[TABLE]

where $v_{\rm rel}$ is the relative velocity between $a$ and $b$.

Appendix B Color Factors

In this appendix, we discuss several processes of colored partner $Y$ and $\bar{Y}$ interaction, and derive the corresponding color factors in the corresponding amplitudes.

We first define the states. For a $Y_{i}\bar{Y}_{j}$ system where $i,j$ are the color indices, the color singlet state is obtained by contraction with $|1\rangle_{ij}=\delta_{ij}/\sqrt{3}$, while the color octet state is obtained by contraction with $|8_{A}\rangle_{ij}=\sqrt{2}(T^{A})_{ij}$. 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 $Y\bar{Y}$ scattering. Clearly, their color indices do not change in such a process, so the color part of the transition operator is simply $\delta_{ik}\delta_{jl}$. The color factor for color singlet to singlet transition is then ${}_{ij}\langle 1|\delta_{ik}\delta_{jl}|1\rangle_{kl}=1$. The color factor for color octet to octet transition is then ${}_{ij}\langle 8_{A}|\delta_{ik}\delta_{jl}|8_{B}\rangle_{kl}=\delta_{AB}$. This is a trivial example.

Second, we consider a gluon emission in $Y\bar{Y}$ 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 $Y$ particle is $(T^{B})_{ik}\delta_{jl}$. And the color factor for octet to singlet transition is ${}_{ij}\langle 1|(T^{B})_{ik}\delta_{jl}|8_{A}\rangle_{kl}=\delta_{AB}/\sqrt{6}$. Squaring the amplitude and sum over $A,B$ yields the color factor 4/3, which is used for deriving Eq. (12).

Third, we consider the Coulomb potential between $Y$ and $\bar{Y}$ generated by a gluon exchange. In this case the operator is $(T^{A})_{ik}(T^{A})^{*}_{jl}=(T^{A})_{ik}(T^{A})_{lj}=\frac{1}{2}\delta_{ij}\delta_{kl}-\frac{1}{6}\delta_{ik}\delta_{jl}$. For the Coulomb potential between $Y\bar{Y}$ in a color singlet, we have the color factor, ${}_{ij}\langle 1|(T^{A})_{ik}(T^{A})_{lj}|1\rangle_{kl}=4/3$. For the Coulomb potential between $Y\bar{Y}$ in a color octet, the color factor is, ${}_{ij}\langle 8_{B}|(T^{A})_{ik}(T^{A})_{lj}|8_{C}\rangle_{kl}=-\delta_{BC}/6$. This explains why $Y\bar{Y}$ in a color singlet (octet) state are attractive (repulsive) to each other, as well as the normalization of the gluon-exchange Coulomb potentials.

Appendix C Capture into Bound States: Gluon Radiated from $Y,\bar{Y}$

In the process $Y_{i}\bar{Y}_{j}\to Y_{i^{\prime}}\bar{Y}_{j^{\prime}}$ ($i,j,i^{\prime},j^{\prime}$ are the color indices), the effective Hamiltonian for an ultra-soft gluon (with color index $C$) radiated from $Y$ or $\bar{Y}$ is

[TABLE]

where $\vec{k}$ is the relative momentum operator between $Y$ and $\bar{Y}$, and $G$ is the gluon field. The matrix element of $H_{int}^{(1)}$ between the final bound state (color singlet) and the initial scattering state (color octet, carrying adjoint color index $D$) is then (see the first two diagrams in Fig 7),

[TABLE]

where $\mu=m_{Y}/2$, and the color factor in the front is obtained based on the discussions in appendix B. $\varepsilon^{A}_{\lambda}$ is the polarization vector of the radiated gluon. The wavefunctions $\Psi_{f}$ and $\Psi_{i}$ are the eigenstates of the Hamiltonians $\hat{H}_{1}=\vec{k}^{2}/(2\mu)-4\alpha_{s}^{(f)}/(3r)$, and $\hat{H}_{8}=\vec{k}^{2}/(2\mu)+\alpha_{s}^{(i)}/(6r)$, respectively, with $\hat{H}_{1}|\Psi_{f}\rangle=E_{f}|\Psi_{f}\rangle$ and $\hat{H}_{8}|\Psi_{i}\rangle=E_{i}|\Psi_{i}\rangle$. Here $\alpha_{S}^{(f)}$ is the strong coupling for the bound state, evaluated at the energy scale equal to the inverse Bohr radius, and $\alpha_{S}^{(i)}$ is the strong coupling for the initial scattering state, evaluated at the energy scale equal to the momentum of the initial $Y$ particle in the center-of-mass frame. We define $\Delta V_{1}=\hat{H}_{8}-\hat{H}_{1}$, which yields Eq. (10).

Using $\vec{k}=-i\mu[\vec{r},H_{8}]=-i\mu\left(\vec{r}H_{8}-H_{1}\vec{r}-\Delta V_{1}(r)\vec{r}\rule{0.0pt}{9.95845pt}\right)$, we get

[TABLE]

Appendix D Capture into Bound States via Non-abelian Gauge Interaction

Because the gluons have their own self-interactions, it is also possible to radiate the gluon from a gluon propagator exchanged between $Y,\ \bar{Y}$ particles, as shown by the third diagram of Fig 7. In non-relativistic 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

[TABLE]

where $f^{ABC}$ is the $SU(3)$ group structure constant. The corresponding covariant derivative on the fundamental $Y$ field is, $D_{\mu}Y_{i}=\partial_{\mu}Y_{i}-ig_{S}G^{A}_{\mu}T^{A}_{ij}Y_{j}$. With these, we calculate the amplitude for the third diagram of Fig 7,

[TABLE]

where $q$ is the momentum running on the gluon propagator.

Going back to the coordinate space, we find this leads to the following interacting Hamiltonian

[TABLE]

Here the strong coupling $\alpha_{S}$ arise from the coupling of $Y$ and $\bar{Y}$ 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 $\alpha_{s}$ at the energy scale which is the average of the inverse Bohr radius and initial $Y$ particle momentum.

Sandwiching $H_{int}^{(2)}$ between the bound and scattering states, we get a new contribution to the transition matrix element

[TABLE]

Summing up Eqs. (30) and (34), we get the total transition matrix element,

[TABLE]

where $\Delta V_{2}(r)=3\alpha_{S}/(2r)$ is also given by Eq. (11). We further derive the capture cross section Eq. (12) based on the above matrix element.

Appendix E Bound States and Projection Operators

In this appendix, we define the bound states made of scalar or fermion color particle and its anti-particle, and introduce the corresponding projection operator for calculating the $S$-wave bound state decay processes.

For a scalar color particle $Y$ (stop or sbottom), we define its plane wave state to be $|\vec{k}\rangle=\sqrt{2\omega_{k}}a^{\dagger}_{\vec{k}}|0\rangle$. Under this convention, $\langle\vec{k}|\vec{k}^{\prime}\rangle=2\omega_{k}(2\pi)^{3}\delta^{3}(\vec{k}-\vec{k}^{\prime})$. Then the bound state at rest made of $Y\bar{Y}$ is defined as

[TABLE]

where $i,j$ are the color indices, and $\tilde{\Psi}(k)$ is the Fourier transformation of the bound state wavefunction in the momentum space, $\tilde{\Psi}(k)=\int d^{3}x\Psi(x)e^{i\vec{k}\cdot\vec{x}}$. It is straightforward to verify that the inner product of two bound state at rest is, $\langle B|B\rangle=4m_{Y}(2\pi)^{3}\delta^{3}(0)$, which agrees with the definition for fundamental particles. When calculating the annihilation decay of a bound state made of $Y\bar{Y}$, we first calculate the amplitude for $Y\bar{Y}\to{\rm final}$, and then use the $Y\bar{Y}$ fields to act on the above bound state. Effectively, the amplitude for bound state decay can be obtained by contracting the ($Y\bar{Y}\to{\rm final}$) amplitude with the following projection operator,

[TABLE]

Similarly, for a bound state made of colored fermion partner and its antiparticle (such as the top partner), the $\eta_{T}$ and $\Upsilon_{T}$ states are defined as

[TABLE]

To calculate the amplitude for bound state decay in this case, we first calculate the amplitude for $Y\bar{Y}\to{\rm final}$ with the external $Y\bar{Y}$ spinors amputated, then contract and trace over it in together with the projector operators of the form

[TABLE]

where $P^{\mu}$ is the four momentum of the decaying bound state.