aainstitutetext: PRISMA Cluster of Excellence & Mainz Institute for
Theoretical Physics, Johannes Gutenberg University, 55099 Mainz, Germany

# A Sommerfeld Toolbox for Colored Dark Sectors

Sonia El Hedri a Anna Kaminska a Maikel de Vries<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present analytical formulae for the Sommerfeld corrections to the
annihilation of massive colored particles into quarks and gluons through the
strong interaction. These corrections are essential to accurately compute the
dark matter relic density for coannihilation with colored partners. Our
formulae allow to compute the Sommerfeld effect, not only for the lowest term
in the angular momentum expansion of the amplitude, but for all orders in the
partial wave expansion. In particular, we carefully account for the effects of
the spin of the annihilating particle on the symmetry of the two-particle wave
function. This work focuses on strongly interacting particles of arbitrary
spin in the triplet, sextet and octet color representations. For typical
velocities during freeze-out, we find that including Sommerfeld corrections on
the next-to-leading order partial wave leads to modifications of up to ten to
twenty percent on the total annihilation cross section. Complementary to QCD,
we generalize our results to particles charged under an arbitrary unbroken
$SU(N)$ gauge group, as encountered in dark glueball models. In connection
with this paper a Mathematica notebook is provided to compute the Sommerfeld
corrections for colored particles up to arbitrary order in the angular
momentum expansion.

††preprint: MITP/16-135††arxiv: 1612.02825

##  1 Introduction

Sommerfeld corrections Sommerfeld:1931aa  through long range interactions play
a critical role in a plethora of thermal dark matter scenarios. Affecting
primarily particles with low velocity, they have been shown, for instance, to
tremendously enhance the dark matter annihilation rate at the galactic center
in various models Hisano:2004ds ; ArkaniHamed:2008qn ; Lattanzi:2008qa . In
particular, the predicted enhancement of the annihilation rate in the galactic
center for pure wino dark matter has allowed to strongly restrict this
supersymmetric scenario Cohen:2013ama ; Fan:2013faa ; Hryczuk:2014hpa ;
Beneke:2016jpw . For multi-TeV weakly interacting dark matter, Sommerfeld
corrections also typically lead to order one modifications of its relic
density Hisano:2006nn ; Hryczuk:2010zi ; Beneke:2014gja ; Beneke:2016ync ,
often significantly weakening the upper bound on the dark matter mass derived
from the Planck measurement Ade:2015xua .

Even in the sub-TeV regime, Sommerfeld corrections become significant in
models involving long range interactions with order one couplings. In
particular, a wide range of dark matter models — such as supersymmetry or
simplified models of coannihilation Baker:2015qna ; Buschmann:2016hkc  —
involve strongly interacting particles in the dark sector. Although the strong
interaction is short ranged at low energies, in the early universe the non-
relativistic QCD potential can be approximated by a Coulomb potential at tree-
level Strumia:2008cf ; deSimone:2014pda . Strongly interacting dark sector
particles would therefore experience sizable long-range interactions through
gluon exchange. These interactions would in turn significantly affect the
annihilation rate into quarks and gluons for masses as low as
$\mathcal{O}(100$ GeV$)$. Computing this rate accurately is crucial in various
scenarios, such as in models where colored particles can survive until short
before BBN, or models where dark matter coannihilates with a colored partner.
In the latter case, the dark matter depletion will in fact be driven by the
annihilation of its coannihilation partner through strong interaction in most
of the parameter space.

Analytical and numerical computations of the Sommerfeld modified annihilation
rate for heavy colored particles have been carried out in various studies
Ellis:2015vaa ; Ibarra:2015nca ; Strumia:2008cf ; deSimone:2014pda ;
Harz:2014gaa . Notably, reference deSimone:2014pda  introduces a general
method to decompose the QCD potential into a sum of Coulomb potentials for
different possible $SU(3)$ representations of the colored dark sector
particle. However, all the existing results only correctly describe
corrections to the $s$-wave cross sections, while higher order effects are
significant. For uncolored particles, the Sommerfeld effect has been computed
beyond the $s$-wave in Cassel:2009wt ; Cassel:2009pu ; Iengo:2009ni . Yet,
these results are consistent only when the annihilation amplitude is dominated
by a single angular momentum component — typically $s$-wave or $p$-wave. In
addition, extending the aforementioned results to colored particles is non-
trivial.

Aside from $SU(3)$, Sommerfeld corrections for dark sector particles charged
under a general $SU(N)$ gauge group have not been considered in the
literature. These non-perturbative effects related to the dark gauge
interaction can significantly modify the annihilation cross-sections of new
charged particles before their freeze-out or, in the case of $SU(N)$ relics,
impact the predicted indirect detection signal. In face of growing interest in
the cosmological role of new gauge groups Gherghetta:2016bcc ; Buen-
Abad:2015ova ; Ko:2016fcd ; Feng:2011ik ; Boddy:2014qxa ; Boddy:2014yra ;
Harigaya:2016nlg ; Soni:2017nlm , the Sommerfeld effect should be derived and
implemented also in the case of non-SM interactions.

In this paper, we present a robust and general framework to analytically
compute the Sommerfeld corrections for the annihilation of dark sector
particles charged either under QCD or $SU(N)$. Instead of considering only the
leading term in the angular momentum expansion of the amplitude, our approach
operates on its complete partial wave expansion into initial states of
definite orbital angular momentum $l$ and spin $s$. Our study focuses on
extensions of the Standard Model with a SM singlet dark matter candidate and
one heavy new particle $\Phi$ which can be either a scalar, a fermion or a
vector. We first consider scenarios where $\Phi$ is a triplet, sextet or octet
of $SU(3)$ and annihilates into quark and gluon pairs. We then generalize
these results to the case where $\Phi$ is charged under either the fundamental
or the adjoint representation of a dark $SU(N)$ gauge group. We discuss direct
applications of these new results in glueball dark matter scenarios. In a
companion paper ElHedri:2017nny , we perform a general study of the relic
density and collider constraints on dark matter models with a colored
coannihilation partner. In these scenarios, the annihilation of $\Phi$ through
strong interactions drives the dark matter depletion and the derived
constraints on the models do not depend on new physics couplings.

The work is organized as follows. In section 2 we discuss the analytic
derivation of Sommerfeld corrections to annihilation processes for arbitrary
partial waves and with any momentum dependence. In section 3 we review
Sommerfeld corrections for QCD potentials in a manner that is applicable to
annihilation of particles with arbitrary color representation. The approaches
in sections 2 and 3 are orthogonal and can be combined into a general
prescription for the annihilation of colored particles. In section 4 we show
that these Sommerfeld effects are significant for colored dark sectors. In
addition to QCD we discuss the Sommerfeld correction for dark sectors charged
under $SU(N)$ in section 5. We conclude in section 6 and discuss more exotic
colored dark sectors in appendix B.

##  2 Sommerfeld corrections for partial waves

Accurately computing the Sommerfeld corrections for an arbitrary process can
prove a daunting task that often has to be performed numerically.
Annihilations in the dark sector, however, involve heavy particles and can
therefore be studied in the non-relativistic limit. In this limit, the tree-
level amplitude for a given process can be reliably approximated by a partial
wave expansion in the orbital angular momentum $l$ and the spin $s$ in either
the initial or final state. Notably, for a $2\rightarrow 2$ process with two
scalar fields in the initial state, this expansion would be of the form

$\mathcal{M}(p,\theta,\phi)=\sum_{l,m}F_{lm}(p)Y_{lm}(\theta,\phi),$ (1)

where $p$ is the magnitude of the initial state momentum in the center-of-mass
frame, $(\theta,\phi)$ are the scattering angles and $Y_{lm}(\theta,\phi)$ are
the spherical harmonics. Without loss of generality, the radial part of the
amplitude can be expanded in powers of $p$ such that the lowest-order
contribution for a given $l$ is $p^{l}$

$F_{lm}(p)=\sum_{n\geq 0}\alpha_{lmn}\,p^{l+2n}.$ (2)

For amplitudes that are dominated by a single partial wave process, the
Sommerfeld corrections can be expressed as an overall multiplicative factor to
the tree-level cross section

$\sigma_{\mathrm{Sommerfeld}}=S\sigma_{\mathrm{perturbative}}.$ (3)

The rescaling factor $S$ encodes the modification of the transition amplitude
by a distorting potential $V$ (modeling the long-range interactions in the
non-relativistic limit) acting on the initial particle wave functions. For a
Coulomb potential $V=-A/r$ in particular, this factor has a simple analytic
form in the $s$-wave Sommerfeld:1931aa

$S(x)=\frac{\pi x}{1-e^{-\pi x}},\quad x=\frac{A}{\beta}.$ (4)

where $\beta$ is the velocity of the incoming particles in the center-of-mass
frame. Positive $A$ corresponds to an attractive potential which leads to an
enhancement of the perturbative result, while negative $A$ results in a
depletion of the cross section due to the repulsive interaction. Analytical
formulae for the Sommerfeld correction factors for higher waves have been
computed in Cassel:2009wt ; Iengo:2009ni  assuming the amplitude is
proportional to $p^{l}$ for the $l$th partial wave. This has been extended
upon slightly in Cassel:2009pu  allowing for a single term with a momentum
dependence of $p^{l+2n}$ with $n\geq 0$. Here, we extend these results to a
full expansion of the annihilation amplitude into orbital angular momentum and
spin states $(l,s)$ up to an arbitrary $l_{\text{max}}$. In particular, we
allow the different terms of the expansion to coexist and we take higher order
terms in equation (2) into account. In the rest of this section we consider a
Coulomb potential and do not make assumptions about the spin of the initial
state particles.

###  2.1 Partial wave expansion

For a given field $\Phi$, the
$\Phi\,\overline{\Phi}\to\textrm{SM}\,\textrm{SM}$ amplitude can be expanded
into orbital angular momentum and spin states $(l,s)$. The reasons for doing
this expansion are manifold. First, as argued at the beginning of this
section, this expansion can be interpreted as a velocity expansion, which
would provide an accurate approximation of the annihilation amplitude for non-
relativistic particles. Moreover, as we will explain in section 2.2, obtaining
Sommerfeld corrections involves computing the non-relativistic wave function
for the two $\Phi$ scattering states. In our case, this wave function is a
solution of the Schrödinger equation for a Coulomb potential. As shown in
Cassel:2009pu , expanding both the scattering state wave function and the
annihilation amplitude considerably simplifies calculations. This leads to a
set of independent equations for each partial wave and allows to obtain
analytical formulae for the Sommerfeld corrected matrix element
$\mathcal{M}^{(S)}_{ls}$. Note that, since the different $(l,s)$ states are
orthogonal to each other, the final cross section will be of the form

$\sigma^{(S)}\propto\sum_{l,s}\left|\mathcal{M}^{(S)}_{ls}\right|^{2}.$ (5)

Another notable advantage of using a $(l,s)$ decomposition is that for
identical particles the overall form of a given $(l,s)$ state is strongly
constrained by CP conservation. For particles carrying no other quantum
numbers than the ones associated to the Lorentz group, a CP transformation
multiplies the initial or final state wave function by $(-1)^{l+s}$. Only
states with even $l+s$ would therefore have a non-zero amplitude. For colored
particles, on the other hand, the color factor in the amplitude can be
decomposed into two parts, respectively symmetric and antisymmetric under
particle exchange. States with even $l+s$ will be proportional to the
symmetric part while states with odd $l+s$ will be proportional to the
antisymmetric part. For $\Phi\,\overline{\Phi}\to g^{a}\,g^{b}$ in particular,
since the gluons are identical particles, the contributions from states with
even $l+s$ will be proportional to
$\left\\{T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}\right\\}$
while the ones for states with odd $l+s$ will be proportional to
$\left[T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}\right]$ where
$T^{a}_{\mathrm{\textbf{R}}}$ is the generator for the representation R of
$\Phi$. This color factor dependence will allow us to introduce a generic
procedure to decompose the amplitude into definite color states as we will
describe in section 3. The same arguments apply to the case of $SU(N)$, which
will be discussed in section 5.

In what follows, we consider a $2\rightarrow 2$ annihilation process in the
center-of-mass frame and in the spin basis. Without loss of generality we
choose the final state particles to be along the $\hat{\mathbf{z}}$-axis and
denote the scattering angles in the initial state by $(\theta,\phi)$. With
$\left\\{m_{1},m_{2}\right\\}$ and $\left\\{m_{3},m_{4}\right\\}$ being the
individual spin projections on the $\hat{\mathbf{z}}$-axis in the initial and
final states respectively, the total annihilation amplitude is defined as an
element of the transition matrix $T$

$\mathcal{T}_{fi}\left(p,\theta,\phi\right)=\langle
p_{f};00;m_{3}m_{4}|T|p;\theta\phi;m_{1}m_{2}\rangle.$ (6)

Here, $p$ and $p_{f}$ being the magnitudes of the momenta in the initial and
final states respectively. The information about the total spins $s_{1,2,3,4}$
in the initial and final state is omitted here for compactness of notation.
Further details about the computation of the total amplitude — notably our
definitions for the momenta and the polarization vectors — are provided in
appendix A.

Decomposing the initial state into states of definite orbital angular momentum
$(l,l_{z})$, the amplitude can be rewritten as111The spherical harmonics
$Y_{l}^{m}(\theta,\phi)$ are normalized as  $\int
Y_{l}^{m}(\theta,\phi)Y_{l^{\prime}}^{m^{\prime}}(\theta,\phi)\,\mathrm{d}\Omega=\delta_{ll^{\prime}}\delta_{mm^{\prime}}.$

$\mathcal{T}_{fi}\left(p,\theta,\phi\right)=\sum_{l,l_{z}}\langle
p_{f};00;m_{3}m_{4}|T|p;ll_{z};m_{1}m_{2}\rangle\,Y^{l_{z}}_{l}(\theta,\phi).$
(7)

A given $|p_{i};ll_{z};m_{1}m_{2}\rangle$ state can be decomposed into
$|p;ll_{z};ss_{z}\rangle$ states using Clebsch-Gordan coefficients

$|p;ll_{z};m_{1}m_{2}\rangle=\sum_{s,s_{z}}\langle
s_{1}m_{1}s_{2}m_{2}|ss_{z}\rangle|p;ll_{z};ss_{z}\rangle,$ (8)

where $s_{1},s_{2}$ are the total spins of the incoming particles. The total
amplitude can then be written as

$\mathcal{T}_{fi}\left(p,\theta,\phi\right)=\sum_{l,l_{z}}\sum_{s,s_{z}}Y^{l_{z}}_{l}(\theta,\phi)\langle
s_{1}m_{1}s_{2}m_{2}|ss_{z}\rangle\mathcal{M}(p;ll_{z};ss_{z};m_{3}m_{4}).$
(9)

The matrix element $\mathcal{M}(p;ll_{z};ss_{z};m_{3}m_{4})$ corresponds to
the contribution of a single initial state $|p_{i};ll_{z};ss_{z}\rangle$ to
the total amplitude. Knowing $\mathcal{T}_{fi}$, this matrix element can be
computed using

$\displaystyle\mathcal{M}(p;ll_{z};ss_{z};m_{3}m_{4})$
$\displaystyle\equiv\langle p_{f};00;m_{3}m_{4}|T|p;ll_{z};sm_{s}\rangle$ (10)
$\displaystyle=\sum_{m_{1},m_{2}}\langle
s_{1}m_{1}s_{2}m_{2}|ss_{z}\rangle\int\mathrm{d}\Omega\,Y^{l_{z}\,*}_{l}(\theta,\phi)\mathcal{T}_{fi}\left(p,\theta,\phi\right).$

Since the $(l,s)$ components of the amplitude are orthogonal the total cross
section is of the form

$\sigma=\frac{1}{64\pi^{2}s}\frac{1}{\sqrt{1-\frac{4m_{\Phi}^{2}}{s}}}\dfrac{1}{d_{\Phi}^{2}d_{\mathrm{\textbf{R}}}^{2}}\sum_{m_{3},m_{4}}\sum_{l,l_{z}}\sum_{s,s_{z}}\left|\mathcal{M}(p;ll_{z};ss_{z};m_{3}m_{4})\right|^{2},$
(11)

where $d_{\Phi}$ is the number of degrees of freedom of the field $\Phi$ and
$d_{\mathrm{\textbf{R}}}$ is the dimensionality of the color representation of
$\Phi$. Another factor $\frac{1}{2}$ needs to be included for identical final
state particles like two gluons.

As mentioned at the beginning of this section, the amplitude for a given $l$
can be expanded in powers of the magnitude of the incoming particle momentum
$p=\sqrt{\frac{s}{4}-m_{\Phi}^{2}}$, with the lowest-order contribution for a
given $l$ being $\mathcal{O}(p^{l})$. We can therefore write

$\mathcal{M}(p;ll_{z};ss_{z};m_{3}m_{4})=\sum_{n\geq
0}\alpha^{(m_{3},m_{4})}_{ll_{z}ss_{z},n}\,p^{l+2n}.$ (12)

Since the matrix element is now expanded in the momentum and in $l$ we can
apply the Sommerfeld corrections to each of the terms in equation (12). This
will be derived in the next section and the total Sommerfeld-corrected cross
sections can then be obtained by the use of equation (11).

###  2.2 Sommerfeld corrections

The Sommerfeld effect is a non-perturbative phenomenon caused by the
distortion of the scattering amplitude of two particles through long range
interactions. This distortion occurs primarily at low velocities and therefore
can particularly affect non-relativistic particles such as the ones in the
dark sector.

Although non-perturbative, the Sommerfeld effect can be approximately modeled
by considering the limit of Feynman diagrams with an infinite number of
particle exchanges Cassel:2009pu . These diagrams should in general include
all the possible two-particle irreducible interactions, which would make the
computation of the final amplitude particularly cumbersome. For non-
relativistic particles, however, the final amplitude is dominated by ladder
diagrams with an infinite number of one-particle exchange iterations such as
the one shown in figure 1. For a given $2\rightarrow n$ process with the
Sommerfeld effect occurring in the initial state, the amplitude then verifies
the following recursion relation Visinelli:2010vg

$\mathcal{M}^{(S)}_{\beta\alpha}=\mathcal{M}^{0}_{\beta\alpha}+\int
d\gamma\frac{\mathcal{M}^{(S)}_{\beta\gamma}V_{\gamma\alpha}}{E_{\alpha}-E_{\gamma}+i\epsilon},$
(13)

where $\alpha$ and $\beta$ are the initial and final states respectively and
the integral over $\gamma$ represents the sum over all possible intermediate
states. $\mathcal{M}^{0}_{\beta\alpha}$ is the perturbative scattering
amplitude corresponding to the exchange of one particle and $V_{\gamma\alpha}$
is the non-relativistic interaction potential distorting the initial state
$\alpha$.

The interaction potential $V_{\gamma\alpha}$ can be rewritten as

$V_{\gamma\alpha}=\langle
k;\theta_{k}\phi_{k};m_{1}m_{2}|\hat{V}|p;\theta_{p}\phi_{p};m_{a}m_{b}\rangle,$
(14)

where $\left\\{m_{1},m_{2}\right\\}$ and $\left\\{m_{a},m_{b}\right\\}$ are
the $z$-components of the spins of the $\alpha$ and $\gamma$ states
respectively and $p,k=|\mathbf{p}|,|\mathbf{k}|$ are the magnitudes of the
momenta $\mathbf{p}$ and $\mathbf{k}$ in these states. In the rest of this
work, we will focus on a spin-independent spherically symmetric potential
$V(|\mathbf{p}-\mathbf{k}|)$. We can therefore factor out the spin states,
which gives

$\displaystyle V_{\gamma\alpha}$ $\displaystyle=\langle
m_{1}m_{2}|m_{a}m_{b}\rangle\langle
k;\theta_{k}\phi_{k}|\hat{V}|p;\theta_{p}\phi_{p}\rangle$ (15)
$\displaystyle=\delta_{m_{1}m_{a}}\delta_{m_{2}m_{b}}V(|\mathbf{p}-\mathbf{k}|).$

$\Phi$$g$$\overline{\Phi}$SMSM Figure 1: Sommerfeld ladder diagram for the
annihilation of $\Phi$ into Standard Model particles.

For initial and final state spins
$\textbf{m}_{i}=\left\\{m_{1},m_{2}\right\\}$ and
$\textbf{m}_{f}=\left\\{m_{3},m_{4}\right\\}$, the Sommerfeld-corrected
amplitude can be expressed as

$\mathcal{M}^{(S)}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{p})=\mathcal{M}^{0}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{p})+\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\frac{\mathcal{M}^{(S)}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{k})V(|\mathbf{p}-\mathbf{k}|)}{E_{\alpha}-E_{\gamma}+i\epsilon}.$
(16)

In the non-relativistic limit, the denominator can be rewritten as

$(E_{\alpha}-E_{\gamma})^{-1}\approx\mathcal{E}-\frac{\mathbf{k}^{2}}{2\mu},$
(17)

where $\mathcal{E}$ is the total energy of the system and $\mu$ is its reduced
mass. For two particles of identical mass $m$, $\mu=\frac{m}{2}$. If the
initial states are off-shell, that is
$\frac{\mathbf{\widetilde{p}}^{2}}{2\mu}\neq\mathcal{E}$, we can define
Iengo:2009ni

$\Phi_{\textbf{m}_{f}\textbf{m}_{i}}(\widetilde{\mathbf{p}})=\frac{\mathcal{M}^{(S)}_{\textbf{m}_{f}\textbf{m}_{i}}(\widetilde{\mathbf{p}})}{\frac{\widetilde{\mathbf{p}}^{2}}{2\mu}-\mathcal{E}},$
(18)

which verifies

$\left(\frac{\widetilde{\mathbf{p}}^{2}}{2\mu}-\mathcal{E}\right)\Phi_{\textbf{m}_{f}\textbf{m}_{i}}(\widetilde{\mathbf{p}})=\mathcal{M}^{0}_{\textbf{m}_{f}\textbf{m}_{i}}(\widetilde{\mathbf{p}})-\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\Phi_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{k})V(|\widetilde{\mathbf{p}}-\mathbf{k}|).$
(19)

In position space (we go from $\widetilde{\mathbf{p}}$ to $\mathbf{r}$) this
becomes

$\left(\frac{-\mathbf{\nabla}^{2}}{2\mu}+V(r)-\mathcal{E}\right)\widetilde{\Phi}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{r})=U^{0}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{r}),$
(20)

where

$U^{0}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{r})=\int\frac{d^{3}q}{(2\pi)^{3}}e^{i\mathbf{r}\cdot\mathbf{q}}\mathcal{M}^{0}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{q}).$
(21)

The final amplitude can now be computed by putting the initial states back on-
shell

$\mathcal{M}^{(S)}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{p})=\mathrm{lim}_{\widetilde{\mathbf{p}}\rightarrow\mathbf{p}}\left(\frac{\widetilde{\mathbf{p}}^{2}}{2\mu}-\mathcal{E}\right)\Phi_{\textbf{m}_{f}\textbf{m}_{i}}(\widetilde{\mathbf{p}})\quad\quad\mathrm{with}\quad\frac{\mathbf{p}^{2}}{2\mu}=\mathcal{E},$
(22)

which leads to Iengo:2009ni ; Cassel:2009wt

$\mathcal{M}^{(S)}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{p})=\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\,\mathcal{M}^{0}_{\textbf{m}_{f}\textbf{m}_{i}}(\mathbf{q})\,\phi_{\mathbf{p}}(\mathbf{q}),$
(23)

where $\phi_{\mathbf{p}}(\mathbf{q})$ obeys the traditional Schrödinger
equation in position space

$\left(\frac{-\nabla^{2}}{2\mu}+V(r)-\frac{p^{2}}{2\mu}\right)\widetilde{\phi}_{\mathbf{p}}(\mathbf{r})=0\quad\quad\mathrm{with}\quad\frac{p^{2}}{2\mu}\equiv\mathcal{E}.$
(24)

For a potential of the form $V(|\mathbf{p}-\mathbf{q}|)$, the wave function
can be rewritten as

$\phi_{\mathbf{p}}(\mathbf{q})=\phi(p,q,\hat{\mathbf{p}}\cdot\hat{\mathbf{q}}),$
(25)

and can therefore be expanded in Legendre polynomials and in spherical
harmonics

$\displaystyle\phi_{\mathbf{p}}(\mathbf{q})$
$\displaystyle=\sum_{l}\frac{2l+1}{4\pi}\,F_{l}(p,q)P_{l}(\hat{\mathbf{p}}\cdot\hat{\mathbf{q}})$
(26)
$\displaystyle=\sum_{l,l_{z}}F_{l}(p,q)Y_{l}^{l_{z}\,*}(\theta_{q},\phi_{q})Y_{l}^{l_{z}}(\theta_{p},\phi_{p}).$

As shown in equation (9), the perturbative amplitude can be expanded in
spherical harmonics as well

$\mathcal{M}^{0}_{\textbf{m}_{i}\textbf{m}_{f}}(\mathbf{q})=\sum_{l,l_{z}}\sum_{s,s_{z}}\langle
s_{1}m_{1}s_{2}m_{2}|ss_{z}\rangle\mathcal{M}(q;ll_{z};ss_{z};\textbf{m}_{f})\,Y^{l_{z}}_{l}(\theta_{q},\phi_{q}).$
(27)

Injecting equations (26) and (27) into equation (23), the Sommerfeld corrected
matrix element can then be decomposed as

$\displaystyle\mathcal{M}^{(S)}_{\textbf{m}_{i}\textbf{m}_{f}}(\mathbf{p})$
$\displaystyle=\sum_{l,l_{z}}\sum_{s,s_{z}}\sum_{l^{\prime},l_{z}^{\prime}}\int\frac{\mathrm{d}^{3}q}{(2\pi)^{3}}\,\langle
s_{1}m_{1}s_{2}m_{2}|ss_{z}\rangle\mathcal{M}(q;ll_{z};ss_{z};\textbf{m}_{f})\,F_{l^{\prime}}(p,q)$
(28) $\displaystyle\qquad\qquad\qquad\quad\times
Y^{l_{z}}_{l}(\theta_{q},\phi_{q})Y^{l_{z}^{\prime}}_{l^{\prime}}(\theta_{q},\phi_{q})Y^{l_{z}^{\prime}}_{l^{\prime}}(\theta_{p},\phi_{p})$
$\displaystyle=\sum_{l,l_{z}}\sum_{s,s_{z}}\langle
s_{1}m_{1}s_{2}m_{2}|ss_{z}\rangle\\!\int\\!\frac{q^{2}\mathrm{d}q}{2\pi^{2}}\,\mathcal{M}(q;ll_{z};ss_{z};\textbf{m}_{f})\,F_{l}(p,q)\,Y^{l_{z}}_{l}(\theta_{p},\phi_{p}).$

Here, in the last line we used the orthogonality relations for the spherical
harmonics. The Sommerfeld corrected amplitude for a given $(l,l_{z},s,s_{z})$
state then takes the following simple form

$\mathcal{M}^{(S)}_{ll_{z};ss_{z};\textbf{m}_{f}}(p)=\int\frac{q^{2}\mathrm{d}q}{2\pi^{2}}\,\mathcal{M}(q;ll_{z};ss_{z};\textbf{m}_{f})\,F_{l}(p,q).$
(29)

Using equation (12), we can re-express this amplitude as

$\mathcal{M}^{(S)}_{ll_{z};ss_{z};\textbf{m}_{f}}(p)=\sum_{n\geq
0}\alpha^{\textbf{m}_{f}}_{ll_{z};ss_{z};n}\int\frac{\mathrm{d}q}{2\pi^{2}}\,q^{l+2n+2}\,F_{l}(p,q).$
(30)

As shown in Cassel:2009pu , the integrals can be rewritten as functions of the
derivatives of the radial components of the wave function $R_{pl}(r)$

$\int\frac{\mathrm{d}q}{2\pi^{2}}\,q^{l+2n+2}\,F_{l}(p,q)=\frac{2^{n}\,n!(2l+2n+1)!!}{(-1)^{n}(-i)^{l}(l+2n)!}\,\left.\frac{\partial^{l+2n}R_{pl}(r)}{\partial
r^{l+2n}}\right|_{r=0}.$ (31)

For a Coulomb potential $V=-A/r$, the radial components of the wave function
can be computed analytically and are equal to

$R_{pl}(z;x)=e^{\frac{\pi
x}{4}}e^{\frac{-iz}{2}}z^{l}\sum_{j=0}^{\infty}\frac{\Gamma\left(1+\frac{ix}{2}+l+j\right)}{(2l+1+j)!}\frac{(iz)^{j}}{j!},$
(32)

where $z=2rp$, $x=Am/p$ and with $p$ and $m$ the momentum and mass of the
incoming particles. In our study, since we consider strong interactions, $A$
will be proportional to the QCD coupling $\alpha_{s}$ or the $SU(N)$ coupling
$\alpha_{N}$.

Using the expression given in equation (32), we can then write

$\mathcal{M}^{(S)}_{ll_{z};ss_{z};\textbf{m}_{f}}(p)=\sum_{n\geq
0}\alpha^{\textbf{m}_{f}}_{ll_{z};ss_{z};n}p^{l+2n}\mathcal{C}_{l}(x)\mathcal{D}_{ln}(x).$
(33)

The Sommerfeld factors $\mathcal{C}_{l}(x)$ and $\mathcal{D}_{ln}(x)$ are
given by

$\displaystyle\mathcal{C}_{l}(x)$
$\displaystyle=\frac{1}{(-i)^{l}}\,e^{\frac{\pi
x}{4}}\,\Gamma\left(1+\frac{ix}{2}\right)\prod_{b=1}^{l}\left(1+\frac{ix}{2b}\right)$
(34) $\displaystyle\mathcal{D}_{ln}(x)$
$\displaystyle=\frac{n!(2l+2n+1)!}{(l+n)!}\sum_{j=0}^{2n}\frac{(-2)^{j}(l+j)!}{j!(2n-j)!(2l+j+1)!}\left[\prod_{b=l+1}^{l+j}\left(1+\frac{ix}{2b}\right)\right],$

where $\mathcal{C}_{l}(x)$ is the correction to the amplitude for a
perturbative matrix element of the form
$p^{l}Y_{l}^{l_{z}}(\theta_{p},\phi_{p})$. Note here that
$\mathcal{D}_{l0}(x)=1$ by construction. The Sommerfeld-corrected squared
matrix element for an $(l,l_{z},s,s_{z})$ initial state as given in equation
(11) can then be written as

$\left|\mathcal{M}^{(S)}_{ll_{z};ss_{z};\textbf{m}_{f}}(p)\right|^{2}=S_{l}(x)\sum_{n,n^{\prime}}\alpha^{\textbf{m}_{f}}_{ll_{z};ss_{z};n}\left(\alpha^{\textbf{m}_{f}}_{ll_{z};ss_{z};n^{\prime}}\right)^{*}\mathcal{D}_{ln}(x)\mathcal{D}_{ln^{\prime}}^{*}(x)\,p^{2(l+n+n^{\prime})},$
(35)

where

$S_{l}(x)=|\mathcal{C}_{l}(x)|^{2}=\frac{\pi x}{1-e^{-\pi
x}}\prod_{b=1}^{l}\left(1+\frac{x^{2}}{4b^{2}}\right)$ (36)

is the Sommerfeld correction for a Coulomb potential and for a perturbative
amplitude of the form $p^{l}Y_{l}^{l_{z}}(\theta_{p},\phi_{p})$ Cassel:2009wt
; Iengo:2009ni . Here, we used $\left|\Gamma(1+ib)\right|=\sqrt{\pi
b\,\mathrm{csch}(\pi b)}$. Note that, in equation (35), since higher order
terms are taken into account in the momentum expansion of the perturbative
amplitude, the Sommerfeld corrections can no longer be factored out. The total
Sommerfeld-corrected cross section is then obtained by plugging equation (35)
into equation (11).

###  2.3 Convergence and strategy

The Sommerfeld corrections as given in equation (35) depend on $l$, $n$ and on
inverse powers of the velocity through $x$. In the perturbative regime, the
angular momentum expansion and velocity expansion of the cross section are
closely related. For a given angular momentum $l$, the lowest-order term of
the perturbative amplitude is at best $\mathcal{O}(v^{l})$ or equivalently
$\mathcal{O}(p^{l})$. This relation is however lost when incorporating the
Sommerfeld corrections. As shown in equation (34), at low velocity, the
Sommerfeld factor for a given $(l,n)$ is $\mathcal{O}(p^{-l-2n-\frac{1}{2}})$.
For a momentum expansion of the perturbative amplitude of the form

$\mathcal{M}_{ll_{z};ss_{z};\textbf{m}_{f}}^{0}(p)=\sum_{n\geq
0}\alpha_{ll_{z};ss_{z};n}^{\textbf{m}_{f}}p^{l+2n},$ (37)

the convergence in the momentum is then jeopardized by the Sommerfeld factor.
The lowest-order term of the momentum expansion of the Sommerfeld-corrected
amplitude given in equation (33) becomes

$\displaystyle\mathcal{M}_{ll_{z};ss_{z};\textbf{m}_{f}}^{(S)}(p)$
$\displaystyle=\sqrt{\frac{\pi Am}{p}}\sum_{n\geq
0}(-1)^{l+n}\alpha_{ll_{z};ss_{z};n}^{\textbf{m}_{f}}m^{l+2n}\,\frac{A^{l+2n}}{2^{l}(l+n)!}\,\frac{n!}{(2n)!}+\mathcal{O}(p^{\frac{1}{2}})$
(38) $\displaystyle=\sqrt{\frac{\pi Am}{p}}\sum_{n\geq
0}\widetilde{\alpha}_{ll_{z};ss_{z};n}^{\textbf{m}_{f}}\,\frac{A^{l+2n}}{2^{l}(l+n)!}\,\frac{n!}{(2n)!}+\mathcal{O}(p^{\frac{1}{2}}),$

where
$\widetilde{\alpha}_{ll_{z};ss_{z};n}^{\textbf{m}_{f}}\equiv(-1)^{l+n}\alpha_{ll_{z};ss_{z};n}^{\textbf{m}_{f}}m^{l+2n}$
is dimensionless. For any value of the orbital angular momentum $l$, the
Sommerfeld-corrected amplitude can then contain terms of order
$p^{-\frac{1}{2}}$. The convergence of the $(l,n)$ expansion of the cross
section is now ensured by the factorial and $2^{l}$ terms as well as by the
powers of $A$ since $A<1$. Hence, the convergence is now in the orbital
angular momentum $l$ instead of the velocity. Nonetheless, due to its
factorial nature the convergence of the corrected cross section is at least as
fast as the one of the perturbative cross section with $l$ and $n$. In fact,
this non-trivial result ensures that the application of the Sommerfeld effect
is a self-consistent procedure.

Since the angular momentum and velocity expansions of the Sommerfeld-corrected
cross section are unrelated, we adopt the following strategy when calculating
Sommerfeld corrections:

  1. 1.

Choose a maximal value $l_{\mathrm{max}}$ for the angular momentum expansion
of both the perturbative and Sommerfeld-corrected cross sections. The choice
for $l_{\mathrm{max}}$ determines the degree of precision for both expansions
according to equations (12) and (38).

  2. 2.

For each value of $l$, include all expansion terms from equation (35) with
$n,n^{\prime}$ satisfying $n+n^{\prime}+l\leq l_{\mathrm{max}}$. This way, the
highest order terms in this expansion are always
$\mathcal{O}(p^{2l_{\mathrm{max}}-1})$. This requirement ensures the
consistency of the expansion of the perturbative cross section in powers of
the incoming momentum.

  3. 3.

Finally, the total Sommerfeld-corrected cross section is obtained by injecting
equation (35) into equation (11).

In this procedure the perturbative amplitude is fully expanded up to
$p^{l_{\mathrm{max}}}$ and the perturbative cross section up to
$p^{2l_{\mathrm{max}}-1}$. Applying Sommerfeld corrections to this expansion
gives an angular momentum expansion of the final cross section up to
$l_{\mathrm{max}}$. In sections 3 and 5 we describe how to embed non-Abelian
gauge theories into this formalism. The results of applying this procedure to
the annihilation of colored particles are shown in section 4.2.

##  3 Sommerfeld corrections for QCD

In the previous section we have computed analytic expressions for the
Sommerfeld corrections of processes with arbitrary partial waves and momentum
dependence. This derivation is based on a Coulomb potential, while the
interactions between colored particles are governed by a QCD potential. An
analytic prescription to decompose the QCD potential as a linear combination
of Coulomb potentials has been first described in Strumia:2008cf ;
deSimone:2014pda  for $s$-wave processes. In this section we extend this
derivation to arbitrary partial waves and point out the differences to the
leading order result. This extension allows for a treatment where higher order
partial waves, arbitrary momentum dependence of the amplitude and QCD effects
can all be taken into account. This prescription allows us to derive an
analytic form for the Sommerfeld corrections of the annihilation of colored
states which we apply to the colored dark sector in the next section.

###  3.1 Decomposing the QCD potential

In order to analytically evaluate the Sommerfeld corrections through the
exchange of soft gluons it is necessary to decompose the QCD potential into a
set of Coulomb-like potentials. This is possible due to the fact the higher-
order QCD potential takes the form Fischler:1977yf ; Schroder:1998vy ;
deSimone:2014pda

$V_{\mathrm{QCD}}=C\frac{\alpha_{s}(\hat{\mu})}{r}\left[1+\frac{\alpha_{s}(\hat{\mu})}{4\pi}\left(c_{1}+2c_{2}(\gamma_{E}+\log\hat{\mu}r)\right)\right]\approx
C\frac{\alpha_{s}(\hat{\mu}\approx 1/r)}{r},$ (39)

where $C$ is proportional to the quadratic Casimir. For example for the quark-
antiquark potential $C=\frac{4}{3}$ and the one-loop coefficients are defined
by $c_{1}=\frac{31}{3}-\frac{10}{9}n_{f}$ and $c_{2}=11-\frac{2}{3}n_{f}$,
where $n_{f}$ is the number of active quark flavors at the scale $\hat{\mu}$.
It shows that the QCD potential at higher orders can be approximated as a
simple Coulomb-like form indicated on the right-hand side of equation (39).
Now, as shown in Strumia:2008cf ; deSimone:2014pda , the QCD potential between
two particles of $SU(3)$ representations R and $\mathrm{\textbf{R}}^{\prime}$
can be rewritten as a sum of Coulomb potentials of the form

$V_{\mathrm{\textbf{R}}\otimes\mathrm{\textbf{R}}^{\prime}}=\frac{\alpha_{s}(\hat{\mu})}{r}\sum_{a}T^{a}_{\mathrm{\textbf{R}}}\otimes
T^{a}_{\mathrm{\textbf{R}}^{\prime}}=\frac{\alpha_{s}(\hat{\mu})}{2r}\sum_{\mathrm{\textbf{Q}}}\Big{[}C_{2}(\mathrm{\textbf{Q}})\mathbbm{1}_{\mathrm{\textbf{Q}}}-C_{2}(\mathrm{\textbf{R}})\mathbbm{1}-C_{2}(\mathrm{\textbf{R}}^{\prime})\mathbbm{1}\Big{]},$
(40)

where
$\mathrm{\textbf{R}}\otimes\mathrm{\textbf{R}}^{\prime}=\bigoplus_{\mathrm{\textbf{Q}}}\mathrm{\textbf{Q}}$
and $C_{2}(\mathrm{\textbf{R}})$, $C_{2}(\mathrm{\textbf{R}}^{\prime})$ are
the quadratic Casimir indices for R and $\mathrm{\textbf{R}}^{\prime}$
respectively. Each irreducible Q component of the initial-state wave function
will then evolve independently in its respective potential. It is important to
note here that $\alpha_{s}(\hat{\mu})$ must be evaluated at a much lower scale
than the hard scale of the annihilation process, namely at scales similar to
the momenta of the incoming particles. For clarity reasons we omit the scale
dependence of $\alpha_{s}$ in the rest of this section.

In what follows, we will consider particle-antiparticle annihilation with
$\mathrm{\textbf{R}}=\mathrm{\textbf{3}},\mathrm{\textbf{6}},\mathrm{\textbf{8}}$
and $\mathrm{\textbf{R}}^{\prime}=\overline{\mathrm{\textbf{R}}}$. The
corresponding color decompositions are

$\displaystyle\mathrm{\textbf{3}}\otimes\overline{\mathrm{\textbf{3}}}$
$\displaystyle=\mathrm{\textbf{1}}\oplus\mathrm{\textbf{8}}$ (41)
$\displaystyle\mathrm{\textbf{6}}\otimes\overline{\mathrm{\textbf{6}}}$
$\displaystyle=\mathrm{\textbf{1}}\oplus\mathrm{\textbf{8}}\oplus\mathrm{\textbf{27}}$
$\displaystyle\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}$
$\displaystyle=\mathrm{\textbf{1}}_{\textbf{S}}\oplus\mathrm{\textbf{8}}_{\textbf{A}}\oplus\mathrm{\textbf{8}}_{\textbf{S}}\oplus\mathrm{\textbf{10}}_{\textbf{A}}\oplus\overline{\mathrm{\textbf{10}}}_{\textbf{A}}\oplus\mathrm{\textbf{27}}_{\textbf{S}}.$

The subscripts S and A indicate whether the representation is symmetric or
antisymmetric respectively under the interchange of the two equal
representations R and $\mathrm{\textbf{R}}^{\prime}$. The quadratic Casimir
indices ($C_{2}$) of these representations are given in table (42) along with
the Dynkin indices, defined as
$C(\mathrm{\textbf{R}})\delta^{ab}=\mathrm{tr}\\!\left(T_{\mathrm{\textbf{R}}}^{a}T_{\mathrm{\textbf{R}}}^{b}\right)$.

R | 1 | 3 | 6 | 8 | 10 | 15 | 27 | 64
---|---|---|---|---|---|---|---|---
$C(\mathrm{\textbf{R}})$ | $0$ | $\frac{1}{2}$ | $\frac{5}{2}$ | $3$ | $\frac{15}{2}$ | $10$ | $27$ | $120$
$C_{2}(\mathrm{\textbf{R}})$ | $0$ | $\frac{4}{3}$ | $\frac{10}{3}$ | $3$ | $6$ | $\frac{16}{3}$ | $8$ | $15$
(42)

Injecting equations (41) and (42) into equation (40), we find

$V_{\mathrm{\textbf{3}}\otimes\overline{\mathrm{\textbf{3}}}}=\\!\frac{\alpha_{s}}{r}\\!\left\\{\\!\\!\begin{aligned}
-\frac{4}{3}&\quad(\mathrm{\textbf{1}})\\\
+\frac{1}{6}&\quad(\mathrm{\textbf{8}})\end{aligned}\right.,\quad
V_{\mathrm{\textbf{6}}\otimes\overline{\mathrm{\textbf{6}}}}=\\!\frac{\alpha_{s}}{r}\\!\left\\{\\!\\!\begin{aligned}
-\frac{10}{3}&\quad(\mathrm{\textbf{1}})\\\
-\frac{11}{6}&\quad(\mathrm{\textbf{8}})\\\
+\frac{2}{3}&\quad(\mathrm{\textbf{27}})\end{aligned}\right.,\quad
V_{\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}}=\\!\frac{\alpha_{s}}{r}\\!\left\\{\\!\\!\begin{aligned}
-3&\quad(\mathrm{\textbf{1}}_{\textbf{S}})\\\
-\frac{3}{2}&\quad(\mathrm{\textbf{8}}_{\textbf{A}},\mathrm{\textbf{8}}_{\textbf{S}})\\\
0&\quad(\mathrm{\textbf{10}}_{\textbf{A}},\overline{\mathrm{\textbf{10}}}_{\textbf{A}})\\\
+1&\quad(\mathrm{\textbf{27}}_{\textbf{S}})\end{aligned}\right..$ (43)

For a particle in a color representation
$\mathrm{\textbf{R}}=\mathrm{\textbf{3}},\mathrm{\textbf{6}},\mathrm{\textbf{8}}$,
the particle-antiparticle QCD potential at tree-level can be decomposed into
Coulomb potentials with coupling strengths set by equation (43). The Coulomb
interaction associated to a given irreducible representation Q will affect the
perturbative annihilation process for which the initial state is in the same
color representation. Computing the Sommerfeld effect for a given annihilation
process therefore requires decomposing the perturbative cross section
according to the color representation of the particle-antiparticle initial
state. Each color channel of the cross section will then be corrected
independently by its own Coulomb potential order by order in the $(l,s)$
expansion. To obtain the full Sommerfeld-corrected amplitude one has to find
the irreducible representations Q contributing at each partial wave order and
the weight of their relative contribution to the process.

###  3.2 Decomposing perturbative cross sections

In this section, we consider tree-level annihilation of a particle $\Phi$ into
quarks and gluons through the strong interaction

$\Phi\,\overline{\Phi}\to
q_{i}\,\bar{q}_{j}\qquad\mathrm{and}\qquad\Phi\,\overline{\Phi}\to
g^{a}\,g^{b}.$ (44)

Since no new physics couplings are involved, the nature of the diagrams
contributing to the annihilation process only depends on the spin of $\Phi$.
Here, we take $\Phi$ to be either a scalar, a fermion or a vector. The Feynman
diagrams for the different annihilation processes are shown in figure 2. Note
that the prescription in this section and the previous section for decomposing
the QCD potential and cross section is also applicable to other processes. A
few more exotic examples are discussed in appendix B.2 and B.3.

$\Phi$$\overline{\Phi}$$q$$\bar{q}$

$\Phi$$\overline{\Phi}$$g$$g$

$\Phi$$\overline{\Phi}$$g$$g$

$\Phi$$\overline{\Phi}$$g$$g$

$\Phi$$\overline{\Phi}$$g$$g$

Figure 2: Feynman diagrams for the annihilation of $\Phi$ into either a quark
anti-quark pair or a pair of gluons. The annihilating field $\Phi$ may be
scalar, fermion or vector, however, in the case of the fermion the four-point
interaction is absent.

First we discuss the color structure of the amplitude for the annihilation
into a quark anti-quark pair. As shown in figure 2, this process occurs
through a single $s$-channel gluon exchange diagram. The corresponding
amplitude is therefore proportional to the generator for the $SU(3)$
representation R of $\Phi$:

$\mathcal{A}^{a}\big{|}^{i}_{j}\propto(T_{\mathrm{\textbf{R}}}^{a})^{i}_{j},$
(45)

where $a$ is the index of the $s$-channel gluon and the indices $i$ and $j$
run from $1$ to the dimensionality of the R representation,
$d_{\mathrm{\textbf{R}}}$. Only the color octet configuration of the initial
state, matching the representation of the exchanged gluon, will therefore
contribute to the $\Phi\,\overline{\Phi}\rightarrow q_{i}\,\bar{q}_{j}$ cross
section

$\sum_{\mathrm{color}}\left|A_{\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}}\right|^{2}=\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2},$
(46)

where $\sum_{\mathrm{color}}$ runs over all the color indices of the external
particles in the amplitude.222Note that this result only holds when
annihilation occurs through an $s$-channel gluon. In some non-minimal models,
$\Phi$ can also annihilate into a quark anti-quark pair through the
$t$-channel exchange of a new particle. We discuss this scenario in more
detail in appendix B.4. The octet representation $[\mathrm{\textbf{8}}]$ is
antisymmetric for the decomposition of self-conjugate representations, like
the $\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}$, because of the CP nature
of the exchanged gluon Kauth:2011vg .

Annihilation processes into gluons have a more complex color structure. As can
be seen in figure 2, four different processes now contribute to the
annihilation cross section, each with a different kinematics. The amplitudes
for all of these processes, however, will be proportional to a linear
combination of $T^{a}_{\mathrm{\textbf{R}}}\,T^{b}_{\mathrm{\textbf{R}}}$ and
$T^{b}_{\mathrm{\textbf{R}}}\,T^{a}_{\mathrm{\textbf{R}}}$ where $a,b$ are the
color indices of the final state gluons. In full generality, the amplitude can
then be written as

$\mathcal{A}^{ab}\big{|}^{i}_{j}=\alpha\left\\{T_{\mathrm{\textbf{R}}}^{a},T_{\mathrm{\textbf{R}}}^{b}\right\\}^{i}_{j}+\beta\left[T_{\mathrm{\textbf{R}}}^{a},T_{\mathrm{\textbf{R}}}^{b}\right]^{i}_{j},$
(47)

where $\alpha,\beta$ are factors that contain the kinematic dependence. As
underlined in section 2.1, this expression drastically simplifies when the
amplitude is expanded into $(l,s)$ states. For a given $(l,s)$ initial state,
CP conservation enforces

$\mathcal{A}^{ab}\big{|}^{i}_{j}=(-1)^{l+s}\mathcal{A}^{ba}\big{|}^{i}_{j}.$
(48)

The annihilation amplitude will therefore be proportional to the
anticommutator of the $T^{a}_{\mathrm{\textbf{R}}}$ for even $l+s$ and to the
commutator for odd $l+s$.333Note that this result differs from the one in
appendix A of deSimone:2014pda  and from Ibarra:2015nca ; Giacchino:2015hvk
which assume proportionality of the amplitude to the anticommutator for all
values of the angular momentum. This simplification allows us to decompose
amplitudes and therefore cross sections into states of definite color
independently of the kinematics of the process.

We now decompose the $\Phi\,\overline{\Phi}\rightarrow g^{a}\,g^{b}$ amplitude
into contributions from initial state configurations with a definite color. As
in section 3.1, we consider particle-antiparticle annihilation with
$\mathrm{\textbf{R}}=\mathrm{\textbf{3}},\mathrm{\textbf{6}},\mathrm{\textbf{8}}$.
For amplitudes proportional to
$[T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}]$, we can write

$\mathcal{A}^{ab}\big{|}^{i}_{j}\propto[T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}]^{i}_{j}=if^{abc}(T^{c}_{\mathrm{\textbf{R}}})^{i}_{j}.$
(49)

As for annihilation into $q\,\bar{q}$, the amplitudes here are proportional to
linear combinations of the generators of the R representation and therefore
receive contributions from color octet configurations only

$\sum_{\mathrm{color}}\left|A_{\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}}\right|^{2}=\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2}.$
(50)

In here the $[\mathrm{\textbf{8}}]$ is in the antisymmetric representation for
the decomposition of self-conjugate representations because of the CP-odd
nature of the amplitude as described in equation (48).

For terms proportional to
$\left\\{T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}\right\\}$,
the amplitude decomposition depends on R. We derive the coefficients
associated to the different color representations of the initial state by
decomposing the amplitude into irreducible tensors as outlined in
Georgi:1999wka . The details of the decomposition of a given tensor for the
processes and representations we are considering are presented in appendix
B.1. Applying the corresponding results to the
$\left\\{T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}\right\\}^{i}_{j}$
tensor allows to express the amplitude as

$\mathcal{A}^{ab}\big{|}^{i}_{j}=\sum_{\mathrm{\textbf{Q}}}\,[\mathrm{\textbf{Q}}]^{ab}\big{|}^{i}_{j}\,,$
(51)

for
$\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}=\bigoplus_{\mathrm{\textbf{Q}}}\mathrm{\textbf{Q}}$
and where $[\mathrm{\textbf{Q}}]^{ab}$ represents the amplitude associated to
an initial state in the color representation Q. Since the contributions from
the different Q initial states are orthogonal, the squared amplitude will be
of the form

$\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}}\big{|}^{2}=\sum_{\mathrm{\textbf{Q}}}\left[\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{Q}}]\big{|}^{2}\right].$
(52)

For the $\Phi\,\overline{\Phi}\rightarrow g^{a}\,g^{b}$ process that we
consider here, using equations (85), (86) and (87) for terms proportional to
$\left\\{T^{a}_{\mathrm{\textbf{R}}},T^{b}_{\mathrm{\textbf{R}}}\right\\}^{i}_{j}$,
we obtain the magnitude of the contributions from the different color states
to the total amplitude. For
$\mathrm{\textbf{R}}=\mathrm{\textbf{3}},\mathrm{\textbf{6}},\mathrm{\textbf{8}}$,
we have

$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{3}}\otimes\bar{\mathrm{\textbf{3}}}}\big{|}^{2}$
$\displaystyle=\frac{7}{2}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}]\big{|}^{2}=\frac{7}{5}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2}$
(53)
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{6}}\otimes\bar{\mathrm{\textbf{6}}}}\big{|}^{2}$
$\displaystyle=\frac{31}{5}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}]\big{|}^{2}=\frac{155}{49}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2}=\frac{155}{81}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{27}}]\big{|}^{2}$
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}}\big{|}^{2}$
$\displaystyle=6\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}_{\mathrm{\textbf{S}}}]\big{|}^{2}=3\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}_{\mathrm{\textbf{S}}}]\big{|}^{2}=2\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{27}}_{\mathrm{\textbf{S}}}]\big{|}^{2}.$

These results for the triplet and the octet agree with the ones obtained for
the $s$-wave in deSimone:2014pda ; Liew:2016hqo . The results for the sextet
and the more exotic decompositions discussed in appendix B.2 and B.3 are novel
and can also be used to extend the scope of the bound state calculations of
Liew:2016hqo  as described in ElHedri:2017nny .

###  3.3 Sommerfeld corrections

Combining the results from sections 3.1 and 3.2, the Sommerfeld-corrected
cross sections for the annihilation of two colored states in the
representations R and $\overline{\mathrm{\textbf{R}}}$ can be decomposed as

$\sigma^{(S)}=\sum_{\mathrm{\textbf{Q}}}\kappa_{\mathrm{\textbf{Q}}}\,\sigma^{(S)}_{\mathrm{C}}\left[\alpha_{\mathrm{\textbf{Q}}}\right],$
(54)

where
$\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}=\bigoplus_{\mathrm{\textbf{Q}}}\mathrm{\textbf{Q}}$.
$\sigma^{(S)}_{C}\left[\alpha_{\mathrm{\textbf{Q}}}\right]$ is the Sommerfeld-
corrected cross section for a Coulomb potential with coupling strength
$A=\alpha_{\mathrm{\textbf{Q}}}$, which can be computed by combining equations
(11) and (35). $\kappa_{\mathrm{\textbf{Q}}}$ is the relative magnitude of the
contribution of the Q initial state to the annihilation amplitude, defined as

$\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{Q}}]\big{|}^{2}=\kappa_{\mathrm{\textbf{Q}}}\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}}\big{|}^{2}.$
(55)

As described in section 3.2, the $\kappa_{\mathrm{\textbf{Q}}}$ weights depend
not only on the color representation of the initial state, but also on its
$(l,s)$ quantum numbers and on the process considered. Notably, for
$\Phi\,\overline{\Phi}\rightarrow g^{a}\,g^{b}$, states with even and odd
$l+s$ are respectively proportional to the anticommutator and the commutator
of the color generators and therefore have different
$\kappa_{\mathrm{\textbf{Q}}}$ factors. In what follows, we will therefore
consider cross sections associated to an individual $(l,s)$ particle-
antiparticle initial state in the
$\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}$ representation.

Reading off $\alpha_{\mathrm{\textbf{Q}}}$ from equation (43) and
$\kappa_{\mathrm{\textbf{Q}}}$ from equations (46), (50) and (53), for
$\mathrm{\textbf{R}}=\mathrm{\textbf{3}},\mathrm{\textbf{6}},\mathrm{\textbf{8}}$,
the Sommerfeld-corrected cross sections are

$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{3}}\otimes\overline{\mathrm{\textbf{3}}}\to
q\,\overline{q}}$
$\displaystyle=\sigma^{(S)}_{C}\left[-\frac{\alpha_{s}}{6}\right]$ (56)
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{3}}\otimes\overline{\mathrm{\textbf{3}}}\to
g\,g}$
$\displaystyle=\begin{cases}\frac{2}{7}\sigma^{(S)}_{C}\left[\frac{4\alpha_{s}}{3}\right]+\frac{5}{7}\sigma^{(S)}_{C}\left[-\frac{\alpha_{s}}{6}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[-\frac{\alpha_{s}}{6}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}$
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{6}}\otimes\overline{\mathrm{\textbf{6}}}\to
q\,\overline{q}}$
$\displaystyle=\sigma^{(S)}_{C}\left[\frac{11\alpha_{s}}{6}\right]$
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{6}}\otimes\overline{\mathrm{\textbf{6}}}\to
g\,g}$
$\displaystyle=\begin{cases}\frac{5}{31}\sigma^{(S)}_{C}\left[\frac{10\alpha_{s}}{3}\right]+\frac{49}{155}\sigma^{(S)}_{C}\left[\frac{11\alpha_{s}}{6}\right]+\frac{81}{155}\sigma^{(S)}_{C}\left[-\frac{2\alpha_{s}}{3}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[\frac{11\alpha_{s}}{6}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}$
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}\to
q\,\overline{q}}$
$\displaystyle=\sigma^{(S)}_{C}\left[\frac{3\alpha_{s}}{2}\right]$
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}\to
g\,g}$
$\displaystyle=\begin{cases}\frac{1}{6}\sigma^{(S)}_{C}\left[3\alpha_{s}\right]+\frac{1}{3}\sigma^{(S)}_{C}\left[\frac{3\alpha_{s}}{2}\right]+\frac{1}{2}\sigma^{(S)}_{C}\left[-\alpha_{s}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[\frac{3\alpha_{s}}{2}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}.$

The Coulomb cross sections $\sigma^{(S)}_{C}\left[\alpha\right]$ can be
readily obtained by plugging the right value for $\alpha$ into the analytic
expressions in section 2. The final analytic expressions for the Sommerfeld-
corrected cross sections can be found by combining equations (11) and (35).

The results in this section have been based on the assumption that
annihilation always involves initial states of definite color. However, as
argued in Berger:2008ti ; Baer:1998pg , rapid interactions of the annihilating
particles with the gluons in the thermal bath may prevent the initial state to
be in a definite color channel. The importance of this effect is unclear since
the time scale may be of the same order as the Sommerfeld effect. Its impact
on the cross section can be bounded by considering an extreme scenario where
annihilation always involves color-averaged initial states. As mentioned in
section 3.2, for annihilation processes into two quarks or into two gluons
with odd $l+s$, the initial state has to always be a color octet. These
processes are therefore not modified by color-averaging. For annihilation into
gluon pairs with even $l+s$ on the other hand, one has to use the averaged
equivalent of equation (43) for the QCD potential. This new potential can be
straightforwardly obtained from equation (56) by averaging over the different
channels. We then obtain

$V_{\mathrm{\textbf{3}}\otimes\overline{\mathrm{\textbf{3}}}\to
g\,g}^{\mathrm{avg}}=-\frac{11}{42}\frac{\alpha_{s}}{r},\quad\qquad
V_{\mathrm{\textbf{6}}\otimes\overline{\mathrm{\textbf{6}}}\to
g\,g}^{\mathrm{avg}}=-\frac{143}{186}\frac{\alpha_{s}}{r},\quad\qquad
V_{\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}\to
g\,g}^{\mathrm{avg}}=-\frac{1}{2}\frac{\alpha_{s}}{r}.$ (57)

This leads to modified Sommerfeld-correction factors for the annihilation into
two gluons with even $l+s$ as

$\displaystyle\sigma^{(S),\,\mathrm{avg}}_{\mathrm{\textbf{3}}\otimes\overline{\mathrm{\textbf{3}}}\to
g\,g}$ $\displaystyle=\sigma^{(S)}_{C}\left[\frac{11\alpha_{s}}{42}\right]$
(58)
$\displaystyle\sigma^{(S),\,\mathrm{avg}}_{\mathrm{\textbf{6}}\otimes\overline{\mathrm{\textbf{6}}}\to
g\,g}$ $\displaystyle=\sigma^{(S)}_{C}\left[\frac{143\alpha_{s}}{186}\right]$
$\displaystyle\sigma^{(S),\,\mathrm{avg}}_{\mathrm{\textbf{8}}\otimes\mathrm{\textbf{8}}\to
g\,g}$ $\displaystyle=\sigma^{(S)}_{C}\left[\frac{\alpha_{s}}{2}\right].$

In the following section, we assume that the annihilation processes occur
through definite color channels. We emphasize however that the Sommerfeld-
corrected cross sections in the color-averaged scenario can also be readily
calculated using our formalism.

##  4 Annihilation in the colored dark sector

In the previous two sections we described how to analytically calculate
Sommerfeld corrections for the annihilation of colored particles including
higher order partial waves. We are now ready to apply these prescriptions to
actual colored dark sectors. We imagine that the dark sector consists of a
single dark matter particle which is a singlet under the Standard Model gauge
groups. Furthermore the dark sector has a colored particle $\Phi$ with
arbitrary spin — scalar, fermion or vector — and with an arbitrary
representation under $SU(3)$. We then introduce a small coupling between DM
and $\Phi$ ensuring chemical and thermal equilibrium between both particles.
The details and phenomenology of this construction are described in an
upcoming accompanying paper ElHedri:2017nny , here we only focus on the
annihilation of the colored particle $\Phi$. We note that in these types of
constructions the relic abundance is completely determined by the annihilation
rate of the colored particle.

These simple models have been introduced for illustrative purposes. We
emphasize however that the methods detailed in this paper are applicable to
the annihilation of colored particles in any kind of dark sector. In the rest
of this section, we introduce a set of simplified models for $\Phi$ and
compute the associated Sommerfeld corrections.

###  4.1 Simplified models

We consider scenarios where $\Phi$ is either a real or complex scalar, Dirac
or Majorana fermion or a real or complex vector boson. The kinetic and mass
terms for $\Phi=\left\\{S,\psi,V\right\\}$ in the complex scalar, Dirac
fermion and complex vector models are then

$\displaystyle\mathcal{L}_{S}$
$\displaystyle=\left[D_{\mu,ij}^{\vphantom{\mu}}S_{j}\right]^{\dagger}\left[D_{ij}^{\mu}S_{j}\right]-m_{S}^{2}\,S_{i}^{\dagger}S_{i}$
(59) $\displaystyle\mathcal{L}_{\psi}$
$\displaystyle=\bar{\psi}_{i}\not{D}_{ij}\psi_{j}-m_{\psi}\bar{\psi}_{i}\psi_{i}$
$\displaystyle\mathcal{L}_{V}$
$\displaystyle=-\frac{1}{2}{V_{\mu\nu,i}}^{\dagger}V^{\mu\nu}_{i}-ig_{s}{V_{i}^{\mu}}^{\dagger}(T^{a}_{\mathrm{\textbf{R}}})_{ij}V_{j}^{\nu}G_{\mu\nu}^{a}+m_{V}^{2}V_{\mu}^{\dagger}V^{\mu},$

where $i,j$ are color indices and the $T^{a}_{\mathrm{\textbf{R}}}$ matrices
are the generators for the color representation R of $\Phi$. To obtain the
Lagrangians for real scalars, Majorana fermions and real vectors each of the
individual terms need to be multiplied by a factor one half. The covariant
derivatives and field strength are given by

$\displaystyle V^{\mu\nu}_{i}$
$\displaystyle=D^{\mu}_{ij}V^{\nu}_{j}-D^{\nu}_{ij}V^{\mu}_{j}$ (60)
$\displaystyle D_{\mu,ij}$
$\displaystyle=\partial_{\mu}\delta_{ij}-ig_{s}G_{\mu}^{a}(T^{a}_{\mathrm{\textbf{R}}})_{ij}.$

Note that the Lagrangian for vectors can also include anomalous terms
Blumlein:1996qp ; Hewett:1993ks  that we chose not to include in this study.
The implications of using a Stückelberg mass term for vector $\Phi$,
especially on perturbative unitarity, are discussed in our companion paper
ElHedri:2017nny .

We list here the analytic cross sections for the pair-annihilation $\Phi$ to
$q\,\bar{q}$ and $g\,g$. The total annihilation cross sections for
$\Phi=S,\psi,V$ are

$\displaystyle\sigma\\!\left(S\,S\\!\to\\!q\,\bar{q}\right)=$
$\displaystyle\frac{2\pi\alpha_{s}^{2}}{3s}\frac{C_{2}(\mathrm{\textbf{R}})}{d_{\mathrm{\textbf{R}}}}\beta_{S}$
(61) $\displaystyle\sigma\\!\left(S\,S\\!\to\\!g\,g\right)=$
$\displaystyle\frac{2\pi\alpha_{s}^{2}}{3s^{3}}\frac{C_{2}(\mathrm{\textbf{R}})}{d_{\mathrm{\textbf{R}}}\beta_{S}^{2}}\\!\Bigg{[}\\!C_{2}(\mathrm{\textbf{G}})\left(s\beta_{S}(10m_{S}^{2}-s)-24m_{S}^{4}\log\frac{1\\!+\\!\beta_{S}}{1\\!-\\!\beta_{S}}\right)\Bigg{.}$
$\displaystyle\Bigg{.}+6C_{2}(\mathrm{\textbf{R}})\left(s\beta_{S}(s+4m_{S}^{2})+m_{S}^{2}(8m_{S}^{2}-4s)\log\frac{1\\!+\\!\beta_{S}}{1\\!-\\!\beta_{S}}\right)\\!\\!\Bigg{]}$
$\displaystyle\sigma\\!\left(\psi\,\bar{\psi}\\!\to\\!q\,\bar{q}\right)=$
$\displaystyle\frac{2\pi\alpha_{s}^{2}}{3s}\frac{C_{2}(\mathrm{\textbf{R}})}{d_{\mathrm{\textbf{R}}}}\frac{1}{\beta_{\psi}}\left(1+\frac{2m_{\psi}^{2}}{s}\right)$
$\displaystyle\sigma\\!\left(\psi\,\bar{\psi}\\!\to\\!g\,g\right)=$
$\displaystyle-\frac{2\pi\alpha_{s}^{2}}{3s^{3}}\frac{C_{2}(\mathrm{\textbf{R}})}{d_{\mathrm{\textbf{R}}}\beta_{\psi}^{2}}\\!\Bigg{[}\\!C_{2}(\mathrm{\textbf{G}})\left(s\beta_{\psi}(s+5m_{\psi}^{2})-12m_{\psi}^{4}\log\frac{1\\!+\\!\beta_{\psi}}{1\\!-\\!\beta_{\psi}}\right)\Bigg{.}$
$\displaystyle\Bigg{.}+3C_{2}(\mathrm{\textbf{R}})\left(s\beta_{\psi}(s+4m_{\psi}^{2})+(8m_{\psi}^{4}-4m_{\psi}^{2}s-s^{2})\log\frac{1\\!+\\!\beta_{\psi}}{1\\!-\\!\beta_{\psi}}\right)\\!\\!\Bigg{]}$
$\displaystyle\sigma\\!\left(V\,V\\!\to\\!q\,\bar{q}\right)=$
$\displaystyle\frac{\pi\alpha_{s}^{2}}{54s}\frac{C_{2}(\mathrm{\textbf{R}})}{d_{\mathrm{\textbf{R}}}}\beta_{V}\frac{12m_{V}^{4}+20m_{V}^{2}s+s^{2}}{m_{V}^{4}}$
$\displaystyle\sigma\\!\left(V\,V\\!\to\\!g\,g\right)=$
$\displaystyle\frac{2\pi\alpha_{s}^{2}}{9m_{V}^{2}s^{3}}\frac{C_{2}(\mathrm{\textbf{R}})}{d_{\mathrm{\textbf{R}}}\beta_{V}^{2}}\\!\Bigg{[}\\!C_{2}(\mathrm{\textbf{G}})m_{V}^{2}\\!\\!\left(\\!s\beta_{V}\\!\left(10m_{V}^{2}\\!+\\!7s\right)\\!-\\!8\\!\left(3m_{V}^{4}\\!+\\!s^{2}\right)\log\frac{1\\!+\\!\beta_{V}}{1\\!-\\!\beta_{V}}\\!\right)\\!\Bigg{.}$
$\displaystyle\Bigg{.}\\!+\\!2C_{2}(\mathrm{\textbf{R}})\\!\left(\\!s\beta_{V}\\!\left(12m_{V}^{4}\\!+\\!3m_{V}^{2}s\\!+\\!4s^{2}\right)\\!+\\!12\left(2m_{V}^{6}\\!-\\!m_{V}^{4}s\right)\log\frac{1\\!+\\!\beta_{V}}{1\\!-\\!\beta_{V}}\\!\right)\\!\\!\Bigg{]}.$

In these expressions, the phase space factor is defined by
$\beta_{\Phi}=\sqrt{1-\frac{4m_{\Phi}^{2}}{s}}$ and
$C_{2}(\mathrm{\textbf{G}})=N$ is the quadratic Casimir of $SU(N)$. The
annihilation cross sections are the same for real scalars, Majorana fermions
and real vectors. Note that, since we directly introduced a squared mass term
for $\Phi$ in the Lagrangian, the $V\,V\rightarrow q\,\bar{q}$ cross section
grows as $\mathcal{O}(s)$ at large center-of-mass energies. This non-physical
behavior can be corrected by introducing a Higgs-type particle. We discuss the
associated effects on the phenomenology in our companion paper ElHedri:2017nny
.

###  4.2 Sommerfeld-corrected annihilation

This section shows the Sommerfeld corrections to the annihilation of colored
particles for the non-relativistic velocities typical to most thermal dark
matter models. Before freeze-out, dark matter and the particles it is in
thermal equilibrium with are forming a thermal bath of relatively low
temperatures compared to their masses. Around freeze-out, when the rate of the
annihilation processes determines the dark matter relic density, the fraction
$x=m/T$ is usually around $25$. This leads to typical velocities around $0.2$
using the Maxwell-Boltzmann distribution. Since, the contributions from larger
velocities are exponentially suppressed, we study the effects of the
Sommerfeld corrections in the thermally relevant range $0<v<0.5$.

We have implemented the procedure detailed in sections 2 and 3 as well as the
perturbative amplitudes for the models described in section 4 in a Mathematica
notebook that is attached to this paper ElHedri:2016pac . This notebook also
provides an interface to micrOMEGAs Belanger:2014vza ; Barducci:2016pcb  for
the calculation of the Sommerfeld-corrected relic abundance in these models.
Furthermore, note that this notebook can also be readily used to compute the
Sommerfeld effect on amplitudes that are not studied here. The conventions and
definitions used to compute the perturbative amplitudes are detailed in
appendix A.

Figure 3: Ratios of the perturbative (solid lines) and Sommerfeld-corrected
cross sections (dashed lines) expanded up to the $s$-wave (blue), $p$-wave
(orange) and $d$-wave (green) over the exact value of the perturbative cross
section. Due to the mass dependence of
$\alpha_{s}^{\mathrm{Sommerfeld}}(\hat{\mu})$, the Sommerfeld-corrected ratios
are shown as a band corresponding to $500~{}\mathrm{GeV}\leq m_{\Phi}\leq
2500~{}\mathrm{GeV}$. For each of the processes we show the results for a
specific color representation, which is denoted by the subscript on the $\Phi$
fields.

In what follows, we consider the ratios of the partial wave expansions of the
perturbative and Sommerfeld-corrected cross sections up to the $d$-wave over
the exact value of the perturbative cross section. For the perturbative cross
sections, we evaluate the strong coupling $\alpha_{s}(\mu)$ at the scale set
by the mass of the annihilating particles. However when taking the ratio of
the cross sections this mass dependence factors out. When computing the
Sommerfeld corrections the coupling
$\alpha_{s}^{\mathrm{Sommerfeld}}(\hat{\mu})$ must be evaluated at a much
lower scale. This is in accordance with the scale of the soft gluons that are
being exchanged. The scale is of the order of the momenta of the incoming
particles that are annihilating and thus depends on the mass of the
annihilating particles and their velocities. Since the scale dependence of
$\alpha_{s}$ is significant for our range of velocities we use the precise
results for the running of the strong coupling obtained in Olive:2016xmw ;
Prosperi:2006hx .

The results for different annihilation processes are shown in figure 3. To
outline the mass dependence of the Sommerfeld-corrected ratios discussed
before, we plot these ratios as a band for $500~{}\mathrm{GeV}\leq
m_{\Phi}\leq 2500~{}\mathrm{GeV}$. We first notice that, as mentioned in
section 2, in spite of the $\mathcal{O}(v^{-1})$ terms present at large $l$
due to Sommerfeld corrections both the perturbative and Sommerfeld-corrected
cross sections converge at similar speeds with $l$. In particular, for all
processes, the $d$-wave perturbative cross section is indistinguishable from
the exact value up to $v\sim 0.5$. Although for colored vectors $l>0$
contributions to the cross sections are negligible, for colored fermions and
scalars, including higher order contributions leads to sizable modifications
of the total cross section for both the perturbative and the Sommerfeld-
corrected case. Notably, for velocities around $0.2$, which is typical for
many thermal dark matter models, adding the $p$-wave contribution can lead to
modifications of $\mathcal{O}(10\%)$ of the Sommerfeld-corrected cross
section. Although in several models these effects can be mitigated by a
cancellation between the $q\,\bar{q}$ and $g\,g$ contributions, our results
highlight the importance of a rigorous computation of Sommerfeld corrections
for more than one partial wave at a time.

As shown in figure 3, the Sommerfeld corrections can enhance the annihilation
cross section of colored particles by up to a factor of two for typical dark
matter velocities. This enhancement plays a crucial role in the phenomenology
of models with a colored dark sector. In an accompanying paper ElHedri:2017nny
, we show how relic density and collider constraints allow to derive model-
independent bounds for scenarios where dark matter coannihilates with a
colored dark partner.

##  5 $SU(N)$ dark sectors

Beyond the minimal models of dark matter explored in section 4, extensions of
the Standard Model involving exotic non-Abelian gauge groups have been
strongly motivated in many BSM theories. In particular, wide classes of models
such as neutral naturalness Gherghetta:2016bcc , hidden valleys
Juknevich:2009ji ; Juknevich:2009gg , dark radiation Buen-Abad:2015ova ;
Ko:2016fcd  and glueball dark matter Soni:2017nlm ; Feng:2011ik ;
Boddy:2014qxa ; Boddy:2014yra ; Harigaya:2016nlg  often involve dark sector
particles charged under a new $SU(N)$ gauge group. When this $SU(N)$ group is
unbroken, dark gluon exchange between the dark sector particles leads to a
long-range interaction through the same mechanism as the one described in
sections 2 and 3 for colored particles. For sizable values of the dark
$\alpha_{N}$ gauge coupling, this long-range interaction leads to significant
Sommerfeld corrections that can be analytically approximated as in the QCD
scenario. Computing the Sommerfeld effect is especially crucial when
considering classes of models where $SU(N)$ is confining in the present
universe Soni:2017nlm ; Feng:2011ik ; Boddy:2014qxa ; Boddy:2014yra ;
Harigaya:2016nlg . Since in these models particles charged under $SU(N)$ are
responsible for dark matter depletion, the Sommerfeld corrections are expected
to significantly change the dark matter relic abundance.

In this section, we extend the methodology outlined in sections 2 and 3 for
QCD to general $SU(N)$ dark sectors. We put special emphasis on the
annihilation of messenger particles charged under both the SM and a dark gauge
group, encountered in large categories of models. We discuss how to combine
the Sommerfeld corrections from both potentials in these scenarios. To
illustrate the relevance of our approach, we compute the Sommerfeld
corrections for the model studied in Soni:2017nlm  that involves dark fermions
charged under both $SU(3)$ and $SU(N)$.

###  5.1 Color decomposition

In this section, we generalize the results derived in section 3 to particles
charged under a new dark gauge group $SU(N)$, either in the fundamental F or
in the adjoint A representation. As before, we consider the self-annihilation
of a particle $\Phi$ into two fermions in the fundamental representation of
$SU(N)$ or into two dark gauge bosons in the adjoint representation of $SU(N)$

$\Phi\,\overline{\Phi}\to
Q_{i}\,\bar{Q}_{j}\qquad\mathrm{and}\qquad\Phi\,\overline{\Phi}\to
G^{a}\,G^{b}.$ (62)

Here, we consider both $Q$ and $G$ to be massless. The procedure for computing
the Sommerfeld corrections for this annihilation process is the same as the
one described for $SU(3)$ in section 3. In particular, the leading order term
of the $SU(N)$ potential can be described by a Coulomb potential obeying
equation (40) and the symmetry constraints on the different partial waves
arising from CP conservation are independent on $N$.

A generalized version of equation (41) can be derived by decomposing the
$\mathrm{\textbf{F}}\otimes\overline{\mathrm{\textbf{F}}}$ and the
$\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}$ products to obtain the
following possible representations for the $\Phi\,\overline{\Phi}$ initial
state

$\displaystyle\mathrm{\textbf{F}}\otimes\overline{\mathrm{\textbf{F}}}$
$\displaystyle=\mathrm{\textbf{1}}\oplus\mathrm{\textbf{A}}$ (63)
$\displaystyle\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}$
$\displaystyle=\mathrm{\textbf{1}}_{\textbf{S}}\oplus\mathrm{\textbf{A}}_{\textbf{A}}\oplus\mathrm{\textbf{A}}_{\textbf{S}}\oplus\mathrm{\textbf{B}}_{\textbf{S}}\oplus\mathrm{\textbf{C}}_{\textbf{A}}\oplus\overline{\mathrm{\textbf{C}}}_{\textbf{A}}\oplus\mathrm{\textbf{D}}_{\textbf{S}}\,.$

A notable difference from the $SU(3)$ case here is the appearance of the
$\mathrm{\textbf{B}}_{\textbf{S}}$ representation for $N\geq 4$. The
representations in equation (63) are associated with the following Young
tableaux

$\ytableausetup{mathmode,boxsize=1em}\mathrm{\textbf{F}}=\raisebox{-1.99168pt}{\ytableau
{} }\qquad\mathrm{\textbf{A}}=\raisebox{14.22636pt}{\ytableau {} & {} \\\ {}
\\\ {} \\\ \none[\raisebox{-1.70717pt}{\large\vdots}] \\\ {}
}\qquad\mathrm{\textbf{B}}=\raisebox{14.22636pt}{\ytableau {} & {} \\\ {} {}
\\\ \none[\raisebox{-1.70717pt}{\large\vdots}] \\\ {}
}\qquad\mathrm{\textbf{C}}=\raisebox{14.22636pt}{\ytableau {} & {} {} \\\ {}
\\\ \none[\raisebox{-1.70717pt}{\large\vdots}] \\\ {}
}\qquad\mathrm{\textbf{D}}=\raisebox{14.22636pt}{\ytableau {} & {} {} {} \\\
{} {} \\\ {} {} \\\ \none[\raisebox{-1.70717pt}{\large\vdots}]
\none[\raisebox{-1.70717pt}{\large\vdots}] \\\ {} {} }\;,$ (64)

where A and D have $N-1$ vertical boxes and B and C have one box less. These
Young tableaux highlight the symmetry properties of the tensors belonging to
the different representations and can therefore be used as guiding tools to
decompose a given amplitude into contributions from different $SU(N)$ initial
states. The dimensionality of all the representations as well as the quadratic
Casimir and Dynkin indices Jeon:2004rk ; Okubo:1977mx  are summarized in the
following table.

R | 1 | F | A | B | C | D
---|---|---|---|---|---|---
$\mathrm{dim}(\mathrm{\textbf{R}})$ | $1$ | $N$ | $N^{2}\\!-\\!1$ | $\frac{1}{4}(N^{4}\\!-\\!2N^{3}\\!-\\!3N^{2})$ | $\frac{1}{4}(N^{4}\\!-\\!5N^{2}\\!+\\!4)$ | $\frac{1}{4}(N^{4}\\!+\\!2N^{3}\\!-\\!3N^{2})$
$C(\mathrm{\textbf{R}})$ | $0$ | $\frac{1}{2}$ | $N$ | $\frac{1}{2}N^{2}(N\\!-\\!3)$ | $\frac{1}{2}N(N^{2}\\!-\\!4)$ | $\frac{1}{2}N^{2}(N\\!+\\!3)$
$C_{2}(\mathrm{\textbf{R}})$ | $0$ | $\frac{N^{2}-1}{2N}$ | $N$ | $2(N\\!-\\!1)$ | $2N$ | $2(N\\!+\\!1)$
(65)

Calling the new gauge coupling $\alpha_{N}$, we can now use equations (40) and
(63) as well as the table in equation (65) to derive the $SU(N)$ Coulomb
potential associated with the different $\Phi\,\overline{\Phi}$
representations

$V_{\mathrm{\textbf{F}}\otimes\overline{\mathrm{\textbf{F}}}}=\\!\frac{\alpha_{N}}{r}\\!\left\\{\\!\\!\begin{aligned}
-\frac{N^{2}-1}{2N}&\quad(\mathrm{\textbf{1}})\\\
\frac{1}{2N}&\quad(\mathrm{\textbf{A}})\end{aligned}\right.,\quad
V_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}}=\\!\frac{\alpha_{N}}{r}\\!\left\\{\\!\\!\begin{aligned}
-N&\quad(\mathrm{\textbf{1}}_{\textbf{S}})\\\
-\frac{N}{2}&\quad(\mathrm{\textbf{A}}_{\textbf{A}},\mathrm{\textbf{A}}_{\textbf{S}})\\\
-1&\quad(\mathrm{\textbf{B}}_{\textbf{S}})\\\
0&\quad(\mathrm{\textbf{C}}_{\textbf{A}},\overline{\mathrm{\textbf{C}}}_{\textbf{A}})\\\
1&\quad(\mathrm{\textbf{D}}_{\textbf{S}})\end{aligned}\right..$ (66)

For the case of $N=3$, this potential reduces to equation (43). For large $N$
the attractive terms increase, whereas the repulsive ones decrease or remain
constant.

The Clebsch-Gordan coefficients for the decomposition of the annihilation
cross sections can now be computed by following exactly the same steps as in
section 3.2. The details of this calculation for the different annihilation
processes as well as for even and odd $l+s$ are given in appendix B.1. As in
the $SU(3)$ case, since the $\Phi\,\overline{\Phi}\to Q_{i}\,\bar{Q}_{j}$
annihilation is mediated by an $s$-channel adjoint gauge boson, only initial
states in the adjoint representation contribute to the total cross section.
For the $\Phi\,\overline{\Phi}\to G^{a}\,G^{b}$ annihilation process, the CP
conservation arguments described in section 3.2 still apply and, as in
equation (50), the squared amplitude for odd $l+s$ can be written as

$\sum_{\mathrm{color}}\left|A_{\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}}\right|^{2}=\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{A}}]\big{|}^{2},$
(67)

for all R. Similarly, for even $l+s$, the decompositions given in equations
(53) for the products of two fundamentals and two adjoints can be generalized
to

$\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{F}}\otimes\overline{\mathrm{\textbf{F}}}}\big{|}^{2}=\frac{N^{2}-2}{2}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}]\big{|}^{2}=\frac{N^{2}-2}{N^{2}-4}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{A}}]\big{|}^{2},$
(68)

and

$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}}\big{|}^{2}$
$\displaystyle=\frac{3}{4}(N^{2}-1)\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}_{\textbf{S}}]\big{|}^{2}$
(69)
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}}\big{|}^{2}$
$\displaystyle=3\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{A}}_{\textbf{S}}]\big{|}^{2}$
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}}\big{|}^{2}$
$\displaystyle=\frac{3(N-1)}{N-3}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{B}}_{\textbf{S}}]\big{|}^{2}$
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}}\big{|}^{2}$
$\displaystyle=\frac{3(N+1)}{N+3}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{D}}_{\textbf{S}}]\big{|}^{2}.$

Note that these results only apply for $N\geq 4$. For $N=3$ the contribution
from the $\mathrm{\textbf{B}}_{\textbf{S}}$ representation goes to zero. In
the case of $N=2$, if $\Phi$ is in the fundamental representation the
$\Phi\,\overline{\Phi}\to G^{a}\,G^{b}$ process occurs only when
$\Phi\,\overline{\Phi}$ is an $SU(2)$ singlet. When $\Phi$ is in the adjoint
representation, only the $\Phi\,\overline{\Phi}$ states in the
$\mathrm{\textbf{1}}_{\textbf{S}}$ and the
$\mathrm{\textbf{D}}_{\textbf{S}}=\mathrm{\textbf{5}}_{\textbf{S}}$
representation will contribute to the $\Phi\,\overline{\Phi}\to G^{a}\,G^{b}$
annihilation cross section. In the large-$N$ limit, on the other hand, we
observe that annihilation to dark gauge bosons occurs dominantly through the
adjoint channel for the annihilation of two fundamentals and splits evenly
into the $\mathrm{\textbf{A}}_{\textbf{S}}$,
$\mathrm{\textbf{B}}_{\textbf{S}}$ and $\mathrm{\textbf{D}}_{\textbf{S}}$
channels for initial state particles in the adjoint representation.

###  5.2 Sommerfeld corrections

We now use the results from section 5.1 as well as the methodology described
in section 3.3 to derive the Sommerfeld correction factors for the
$\Phi\,\overline{\Phi}\to Q_{i}\,\bar{Q}_{j}$ and the
$\Phi\,\overline{\Phi}\to G^{a}\,G^{b}$ annihilation processes. For general
$N$, these factors now read

$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{F}}\otimes\overline{\mathrm{\textbf{F}}}\to
Q\,\overline{Q}}$
$\displaystyle=\sigma^{(S)}_{C}\left[-\frac{\alpha_{N}}{2N}\right]$ (70)
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{F}}\otimes\overline{\mathrm{\textbf{F}}}\to
G\,G}$
$\displaystyle=\begin{cases}\frac{2}{N^{2}-2}\sigma^{(S)}_{C}\left[\frac{(N^{2}-1)\alpha_{N}}{2N}\right]+\frac{N^{2}-4}{N^{2}-2}\sigma^{(S)}_{C}\left[-\frac{\alpha_{N}}{2N}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[-\frac{\alpha_{N}}{2N}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}$
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}\to
Q\,\overline{Q}}$
$\displaystyle=\sigma^{(S)}_{C}\left[\frac{N\alpha_{N}}{2}\right]$
$\displaystyle\sigma^{(S)}_{\mathrm{\textbf{A}}\otimes\mathrm{\textbf{A}}\to
G\,G}$
$\displaystyle=\begin{cases}\frac{4}{3(N^{2}-1)}\sigma^{(S)}_{C}\left[N\alpha_{N}\right]+\frac{1}{3}\sigma^{(S)}_{C}\left[\frac{N\alpha_{N}}{2}\right]&\\\
+\frac{N-3}{3(N-1)}\sigma^{(S)}_{C}\left[\alpha_{N}\right]+\frac{N+3}{3(N+1)}\sigma^{(S)}_{C}\left[-\alpha_{N}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[\frac{N\alpha_{N}}{2}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}.$

The ratios of the $s$, $p$ and $d$-wave annihilation cross sections are shown
in figure 4 for the $\Phi\,\overline{\Phi}\to Q_{i}\,\bar{Q}_{j}$ and
$\Phi\,\overline{\Phi}\to G^{a}\,G^{b}$ processes with $\Phi$ being either a
scalar or a fermion, in either the fundamental or the adjoint representation
of $SU(N)$. As in figure 3, we consider velocity expansions of the cross
section up to the $s$-wave, the $p$-wave and the $d$-wave but this time, we
show the values of these different cross sections for $4\leq N\leq 10$.
Contrary to the $SU(3)$ case, we do not evaluate $\alpha_{N}$ at the scale of
the momenta of the incoming particles and instead set the coupling entering
into the Sommerfeld corrections to be
$\alpha_{N}^{\mathrm{Sommerfeld}}(\hat{\mu})=0.1$. For the typical momenta
considered here, this value is lower than the ones encountered in the QCD
case, thereby leading to conservative estimates of the Sommerfeld effect in
strongly coupled theories.

Figure 4: Ratios of the perturbative (solid lines) and Sommerfeld-corrected
cross sections (dashed lines) expanded up to the $s$-wave (blue), $p$-wave
(orange) and $d$-wave (green) over the total perturbative cross section. We
show the ratios as a band corresponding to $4\leq N\leq 10$ for a specific
process and representation (either F or A), which is denoted by the subscript
on the $\Phi$ fields.

For fundamental particles in the initial state, the Sommerfeld corrections
become negligible in the large $N$ limit. This result can be understood by
noting that, in equation (70), either the effective couplings for the Coulomb
potentials or the coefficients of the $\sigma_{C}$ cross sections are
inversely proportional to powers of $N$. For initial state particles in the
adjoint representation, however, the dominant contributions in the large $N$
limit arise from terms of the form
$\sigma_{C}\left[\frac{N\alpha_{N}}{2}\right]$. In this case, the Sommerfeld
enhancement will therefore grow with $N$ for each partial wave contribution,
as can be observed in figure 4. Note that in this scenario, the Sommerfeld
enhancement is extremely relevant at typical freeze-out velocities and taking
it into account is essential for relic abundance computations.

###  5.3 Messenger particles

One particular scenario often encountered in the literature is the existence
of new particles that are charged both under QCD and under a new $SU(N)$ gauge
group. These particles notably play the roles of messengers between the
Standard Model and the dark sector in hidden valley models Juknevich:2009ji ;
Juknevich:2009gg ; Soni:2017nlm . In this case the non-relativistic potential
for the Sommerfeld effect is the sum of the $SU(3)$ and the $SU(N)$ potential
and Sommerfeld correction factors are modified accordingly Zhang:2013qza . The
total potential is given as

$V=\frac{\alpha_{s}(\hat{\mu})}{r}\sum_{a}T^{a}_{\mathrm{\textbf{R}}}\otimes
T^{a}_{\mathrm{\textbf{R}}^{\prime}}+\frac{\alpha_{N}(\hat{\mu})}{r}\sum_{a}T^{a}_{\mathrm{\textbf{P}}}\otimes
T^{a}_{\mathrm{\textbf{P}}^{\prime}}\,,$ (71)

where R is a representation of $SU(3)$ and P is a representation of $SU(N)$.
Computing this potential for initial states in different $SU(3)\times SU(N)$
representations can be done by applying equation (40) to each of the terms on
the right hand side of equation (71) separately. For a given annihilation
process, the Clebsch-Gordan coefficient for an initial state with given
$SU(3)\times SU(N)$ quantum numbers is the product of the coefficients
corresponding to the $SU(3)$ and the $SU(N)$ representations. These
coefficients can be readily computed using equations (67), (68) and (69). In
what follows, we will apply this procedure to the particular case of particles
charged under the fundamental representations of both $SU(3)$ and $SU(N)$.

###  5.4 Application: bifundamental messengers

In models where the gauge bosons of the dark $SU(N)$ either are dark radiation
or form dark glueballs a connection between the dark sector and the Standard
Model needs to exist to ensure thermal equilibrium. This connection can be
established by introducing messenger particles charged both under QCD and
under the dark $SU(N)$ gauge group. These particles are initially in thermal
equilibrium with the SM and therefore annihilate to SM particles until they
freeze out. When the temperature of the universe drops below the confining
scale of the theory at later times, these messengers form bound states that
decay to dark gauge bosons that ultimately form stable glueball dark matter
candidates. The strength of the messenger annihilation channels to the visible
and dark sectors will therefore set the dark matter relic abundance.

Figure 5: Sommerfeld correction factors for the $s$-wave annihilation cross
sections of fermionic messengers in the fundamental representation of both QCD
and $SU(N)$. The different colors show the relevant annihilation processes and
the different lines represent $N=3$ (solid), $N=5$ (dotted) and $N=10$
(dashed).

In what follows, we consider a fermionic messenger particle $\psi$ charged as
a triplet under QCD and as a fundamental under $SU(N)$. In this scenario,
$\psi$ can annihilate either to $g\,g$, $q\,\bar{q}$, $G\,G$ or $g\,G$, where
$G$ is the massless dark gauge boson for the $SU(N)$ gauge group. The first
two processes occur through the QCD interaction and, since the final states
are $SU(N)$ singlets, the initial $\psi\,\bar{\psi}$ state must also be an
$SU(N)$ singlet. The different QCD representations for $\psi\,\bar{\psi}$ as
well as their corresponding Clebsch-Gordan coefficients are therefore the ones
derived in section 3.3. As outlined in section 5.3, however, the non-
relativistic potential between the two initial state particles will now have
an additional term corresponding to the exchange of dark gluons. Since
$\psi\,\bar{\psi}$ has to be an $SU(N)$ singlet, the new potential will be of
the form

$V=V_{SU(3)}-\frac{N^{2}-1}{2N}\frac{\alpha_{N}}{r},$ (72)

where $V_{SU(3)}$ is given in equation (43). The $\psi\,\bar{\psi}\to G\,G$
process occurs through $SU(N)$ interactions and has been studied in section
5.2. The results from this section can be directly applied to this scenario
with the potential being modified as

$V=-\frac{4}{3}\frac{\alpha_{s}}{r}+V_{SU(N)},$ (73)

where $V_{SU(N)}$ is given in equation (63). Note here that $\psi\,\bar{\psi}$
now is an $SU(3)$ singlet.

Finally, the $\psi\,\bar{\psi}\to g\,G$ process has not been studied before
and has been not been taken into account in Soni:2017nlm . For this
annihilation channel, gauge conservation constrains the $\psi\,\bar{\psi}$
initial state to be in the adjoint representation of both $SU(3)$ and $SU(N)$.
Hence, there is no need to compute any Clebsch-Gordan coefficient and the
potential will now read

$V=\frac{1}{6}\frac{\alpha_{s}}{r}+\frac{1}{2N}\frac{\alpha_{N}}{r}.$ (74)

The Sommerfeld-corrected annihilation cross sections for all these processes
in the $s$-wave can then be expressed as

$\displaystyle\left(\sigma v\right)_{\psi\,\bar{\psi}\to q\,\bar{q}}=$
$\displaystyle\,6\times\frac{\pi\alpha_{s}^{2}}{9Nm_{Q}^{2}}\times
S\left(-\frac{\alpha_{s}}{6\beta}+\frac{N^{2}\\!-\\!1}{2N}\frac{\alpha_{N}}{\beta}\right)$
(75) $\displaystyle\left(\sigma v\right)_{\psi\,\bar{\psi}\to g\,g}=$
$\displaystyle\,\frac{7\pi\alpha_{s}^{2}}{54Nm_{Q}^{2}}\times\left[\frac{2}{7}S\left(\frac{4\alpha_{s}}{3\beta}+\frac{N^{2}\\!-\\!1}{2N}\frac{\alpha_{N}}{\beta}\right)+\frac{5}{7}S\left(-\frac{\alpha_{s}}{6\beta}+\frac{N^{2}\\!-\\!1}{2N}\frac{\alpha_{N}}{\beta}\right)\right]$
$\displaystyle\left(\sigma v\right)_{\psi\,\bar{\psi}\to G\,G}=$
$\displaystyle\,\frac{(N^{2}\\!-\\!1)(N^{2}\\!-\\!2)\pi\alpha_{N}^{2}}{48N^{3}m_{Q}^{2}}$
$\displaystyle\,\times\left[\frac{2}{N^{2}\\!-\\!2}S\left(\frac{4\alpha_{s}}{3\beta}+\frac{N^{2}\\!-\\!1}{2N}\frac{\alpha_{N}}{\beta}\right)+\frac{N^{2}\\!-\\!4}{N^{2}\\!-\\!2}S\left(\frac{4\alpha_{s}}{3\beta}-\frac{1}{2N}\frac{\alpha_{N}}{\beta}\right)\right]$
$\displaystyle\left(\sigma v\right)_{\psi\,\bar{\psi}\to g\,G}=$
$\displaystyle\,\frac{2(N^{2}\\!-\\!1)\pi\alpha_{s}\alpha_{N}}{9N^{2}m_{Q}^{2}}\times
S\left(-\frac{\alpha_{s}}{6\beta}-\frac{1}{2N}\frac{\alpha_{N}}{\beta}\right),$

where the Sommerfeld factor $S(x)$ is given in equation (4). The ratios of
these cross sections over the $s$-wave perturbative cross sections for each
process are shown in figure 5 for $N=3,5,10$ and with
$\alpha_{s}(\hat{\mu})=\alpha_{N}(\hat{\mu})=0.1$. For typical freeze-out
velocities $v\sim 0.2$, the Sommerfeld effect can lead to a factor of two to
eight enhancement of the annihilation cross section for most processes.
Although this enhancement could be slightly mitigated by the reduction of the
cross section for $\psi\,\bar{\psi}\to g\,G$, this reduction is in general
much less pronounced than the enhancement observed for the other processes,
especially as $N$ increases. Taking the Sommerfeld corrections into account
for the annihilation of messenger particles is therefore essential to derive
robust cosmological bounds for the hidden sector models of dark matter
discussed in Soni:2017nlm ; Feng:2011ik ; Boddy:2014qxa ; Boddy:2014yra ;
Harigaya:2016nlg .

##  6 Conclusions

In this work we have derived analytical expressions for the Sommerfeld
corrections of the annihilation of colored particles. These expressions result
from combining two orthogonal procedures: deriving Sommerfeld corrections for
partial waves beyond the leading order and decomposing the QCD potential into
Coulomb potentials. Our results significantly improve on existing literature
and allow to combine higher-order velocity corrections with the QCD nature of
these annihilation processes. These analytical expressions can readily be
applied to any type of annihilation of colored particles in dark sector. The
only necessary step is to expand the annihilation cross sections into states
of definite orbital angular momentum and spin $(l,s)$ and then apply the
correction factors as presented in our work.

For consistently applying Sommerfeld correction factors for higher partial
waves we showed it is necessary to expand the annihilation amplitudes in
$(l,s)$ states. Then one can further expand these states in powers of the
momentum and solve the non-relativistic Schrödinger equation for each of the
states separately. From these solutions one obtains the analytic Sommerfeld-
correction factors for all orders in the partial wave expansion and all powers
of the momentum. We express these results conveniently as the Sommerferld-
correction factor for the $s$-wave times an analytic distortion factor
specific to each term in the partial wave expansion.

The QCD nature of the process poses a challenge for the analytic calculation
of the Sommerfeld corrections which can be overcome by decomposing the
potential into a linear combination of Coulomb potentials. This procedure,
however, depends crucially on the symmetry properties of the amplitude. With
an expansion of the amplitude in $(l,s)$ states these properties become
apparent. The color structure then simplifies and becomes independent of the
kinematics of the process. Then the color-dependent part of the annihilation
amplitude can be treated separately and later combined with the Sommerfeld
corrections for the partial wave components.

Finally, we apply these results to several colored dark sectors with a singlet
dark matter candidate, where the annihilation of the colored states is solely
responsible for setting the relic abundance. We show that for particles of any
spin — scalar, fermion, vector — and in the triplet, sextet or octet
representation of QCD Sommerfeld corrections are sizable. A consistent and
precise inclusion of these effects is therefore essential in understanding the
specific details of a possible colored dark sector. In an accompanying paper
we present the full study of several types of colored dark sectors where we
include precise determination of the relic density and discuss the full
phenomenology of these models.

We also present the first calculation of the Sommerfeld corrections for dark
sectors charged under general $SU(N)$ gauge groups. These corrections are
especially relevant in scenarios where confinement occurs after freeze-out,
since the new gauge group remains unbroken and the gauge coupling is sizable.
Although non-perturbative effects for these models have been previously
overlooked in the literature, we showed that the Sommerfeld corrections can
drastically modify the annihilation cross section of dark sector particles,
and therefore the dark matter relic density. We advocate for taking these
corrections into account in future in-depth studies of these models.

We conclude by emphasizing that the procedure described in this work is not
restricted to the annihilation of identical particles. Notably, our method
also applies to processes like the annihilation of a triplet and an octet of
QCD — for example squark-gluino annihilation in supersymmetry. Henceforth,
Sommerfeld corrections for models with extended dark sectors and multiple
gauge groups can be easily included using the presented formalism.

###### Acknowledgements.

We would like to thank Andrzej Hryczuk, Felix Yu and José Zurita for valuable
discussions. This research is supported by the Cluster of Excellence Precision
Physics, Fundamental Interactions and Structure of Matter (PRISMA-EXC 1098),
by the ERC Advanced Grant EFT4LHC of the European Research Council, and by the
Mainz Institute for Theoretical Physics.

##  Appendix A Partial wave cross sections

This appendix details the conventions used to compute the annihilation cross
sections in the Mathematica notebook attached to this paper ElHedri:2016pac .
In order for the amplitudes to be of the form of equation (6), we work in the
so-called _final frame_ where the momenta of the final state are along the
$z$-axis and the momenta of the initial states are characterized by the angles
$\theta$ and $\phi$. In order for the polarization vectors of the gluons and
of the vector $\Phi$ to be well-defined, we compute the amplitudes in the
helicity basis based on Just:2008xy . In this basis, for the annihilation
processes we consider in the main body of the paper the $\phi$ dependence of
the amplitude is well-known

$\mathcal{A}_{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}}(p,\theta,\phi)\equiv\mathcal{A}_{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}}(p,\theta,\phi=0)e^{i(\lambda_{i}-\lambda_{f})\phi},$
(76)

with $\lambda_{i}=\lambda_{1}-\lambda_{2}$ and
$\lambda_{f}=\lambda_{3}-\lambda_{4}$.444For other processes like the ones
discussed in appendices B.3 and B.4 one needs to take into account the full
$\phi$-dependence of the amplitude. For these processes, we therefore compute
the amplitudes for $\phi=0$ and inject the azimuthal phase factor into the
final expression.

Assuming the quarks to be massless, the momenta for the annihilation of a pair
of $\Phi$’s into quark and gluon pairs are

$\displaystyle p_{1}$
$\displaystyle=(E,p\sin\theta\cos\phi,p\sin\theta\sin\phi,p\cos\theta)$ (77)
$\displaystyle p_{2}$
$\displaystyle=(E,-p\sin\theta\cos\phi,-p\sin\theta\sin\phi,-p\cos\theta)$
$\displaystyle p_{3}$ $\displaystyle=(E,0,0,E)$ $\displaystyle p_{4}$
$\displaystyle=(E,0,0,-E),$

with $E=\sqrt{p^{2}+m^{2}}$. The spinors for a particle of helicity
$\pm\frac{1}{2}$ and mass $m$ moving in the direction $(\theta,\phi)$ are

$\displaystyle u_{+}(p,m,\theta,\phi)$
$\displaystyle=\mathcal{R}(\theta,\phi)\cdot\begin{pmatrix}\sqrt{E-p}\\\ 0\\\
\sqrt{E+p}\\\ 0\end{pmatrix}\quad\quad
u_{-}(p,m,\theta,\phi)=\mathcal{R}(\theta,\phi)\cdot\begin{pmatrix}0\\\
\sqrt{E+p}\\\ 0\\\ \sqrt{E-p}\end{pmatrix}$ (78) $\displaystyle
v_{+}(p,m,\theta,\phi)$
$\displaystyle=\mathcal{R}(\theta,\phi)\cdot\begin{pmatrix}0\\\ -\sqrt{E+p}\\\
0\\\ \sqrt{E-p}\end{pmatrix}\quad
v_{-}(p,m,\theta,\phi)=\mathcal{R}(\theta,\phi)\cdot\begin{pmatrix}\sqrt{E-p}\\\
0\\\ -\sqrt{E+p}\\\ 0\end{pmatrix},$

with

$\mathcal{R}(\theta,\phi)=\begin{pmatrix}\cos\frac{\theta}{2}&-\sin\frac{\theta}{2}e^{-i\phi}&0&0\\\
\sin\frac{\theta}{2}e^{i\phi}&\cos\frac{\theta}{2}&0&0\\\
0&0&\cos\frac{\theta}{2}&-\sin\frac{\theta}{2}e^{-i\phi}\\\
0&0&\sin\frac{\theta}{2}e^{i\phi}&\cos\frac{\theta}{2}\end{pmatrix}.$ (79)

The spinors for a particle moving in the opposite direction are obtained in
Just:2008xy  as well. They are very similar to the expressions in equation
(78) and are given in the attached Mathematica notebook ElHedri:2016pac . When
computing amplitudes involving fermion currents, we define the gamma matrices
in the Weyl basis.

The transverse polarization vectors corresponding to a final state gluon of
momentum $p_{3}$ or to a vector $\Phi$ of momentum $p_{1}$ are

$\epsilon_{\pm}^{(1)}(\theta)=\frac{1}{\sqrt{2}}\left(0,\mp\cos\theta,-i,\pm\sin\theta\right),$
(80)

while the longitudinal polarization vector corresponding to a vector $\Phi$ of
momentum $p_{1}$ is

$\epsilon_{0}^{(1)}(p,m,\theta)=\left(\frac{p}{m},\frac{E}{m}\sin\theta,0,\frac{E}{m}\cos\theta\right).$
(81)

Similarly, the transverse polarization vectors corresponding to a final state
gluon of momentum $p_{4}$ or to a vector $\Phi$ of momentum $p_{2}$ are

$\epsilon_{\pm}^{(2)}(\theta)=\frac{1}{\sqrt{2}}\left(0,\pm\cos\theta,-i,\mp\sin\theta\right),$
(82)

while the longitudinal polarization vector corresponding to a vector $\Phi$ of
momentum $p_{2}$ is

$\epsilon_{0}^{(2)}(p,m,\theta)=\left(-\frac{p}{m},\frac{E}{m}\sin\theta,0,\frac{E}{m}\cos\theta\right).$
(83)

The full $\phi$-dependence for polarization vectors can be found in the
attached Mathematica notebook ElHedri:2016pac . We use Chung:1971ri  to
convert the helicity amplitudes into amplitudes for definite spin states
$m_{1},m_{2},m_{3},m_{4}$ (with $m_{i}$ being the $z$-component of the spin of
particle $i$) using

$\mathcal{A}_{m_{1},m_{2},m_{3},m_{4}}(p,\theta,\phi)\\!=\\!\\!\sum_{\lambda_{1},\lambda_{2}}\\!D^{s_{1}\,*}_{m_{1},\lambda_{1}}\\!(\phi,\theta,-\phi)D^{s_{2}\,*}_{m_{2},-\lambda_{2}}\\!(\phi,\theta,-\phi)\mathcal{A}_{\lambda_{1},\lambda_{2},m_{3},-m_{4}}(p,\theta,\phi)$
(84)

where the $D^{s}_{m,\lambda}$ are the Wigner D-functions. In this formula, we
used $\lambda_{3}=m_{3}$ and $\lambda_{4}=-m_{4}$ for final states, whose
momenta are along the $z$-axis.

##  Appendix B Color decomposition

In this appendix we describe the decomposition of the group structure of the
amplitudes discussed in section 3.2 for QCD and in section 5.1 for $SU(N)$.
Moreover, later in this appendix we discuss the decomposition and Sommerfeld
corrections for more exotic particles present in colored dark sectors.

###  B.1 Amplitude tensor decomposition

In this section we describe how to decompose a colored amplitude into several
channels of definite color. Then using these expressions and using a specific
form for the color part of the amplitude as obtained in section 3.2 we square
the amplitude and find the decomposition of the total cross section. In
principle one can decompose amplitudes which may be any product of
representations of $SU(N)$, however, here we restrict ourselves to
$\mathrm{\textbf{R}}\otimes\overline{\mathrm{\textbf{R}}}$ with
$\mathrm{\textbf{R}}=\mathrm{\textbf{F}}$ (fundamental), S
($\frac{N(N+1)}{2}$-dimensional symmetric) and A (adjoint). The decomposition
for these channels is given in equation (41) for QCD and in equation (63) for
$SU(N)$. More exotic combinations are discussed for QCD in the next two
sections of this appendix. To decompose the amplitudes we base ourself on the
method of tensor decomposition as described in Georgi:1999wka  and use
fundamental indices for all representations. To switch between
$(T_{\mathrm{\textbf{R}}}^{a})^{i}_{j}$ where $i,j$ run from $1$ to
$d_{\mathrm{\textbf{R}}}$ and the fundamental indices one can use the Clebsch-
Gordan coefficients of the representation with respect to fundamentals of
$SU(N)$. We focus purely on the color of $\Phi$ in $A$ and drop the color
dependence of the remaining part in the amplitudes. This method has been put
forward already in deSimone:2014pda  for the triplet and the octet in QCD and
we extend these results to arbitrary $N$. Parts of these calculation have been
done using LieArt Feger:2012bs  and ColorMath Sjodahl:2012nk . For the product
of two fundamentals we can write the tensor product as $A^{k}_{i}=v^{i}w_{j}$
which contains the full color dependence of the total amplitude. We split up
this part of the amplitude as

$\displaystyle A^{i}_{j}$
$\displaystyle=[\mathrm{\textbf{1}}]^{i}_{j}+[\mathrm{\textbf{A}}]^{i}_{j}$
(85) $\displaystyle[\mathrm{\textbf{1}}]^{i}_{j}$
$\displaystyle=\frac{1}{N}\delta^{i}_{j}A^{m}_{m}$
$\displaystyle[\mathrm{\textbf{A}}]^{i}_{j}$
$\displaystyle=A^{i}_{j}-\frac{1}{N}\delta^{i}_{j}A^{m}_{m}.$

In here the indices $i,j=1,\cdots,N$ represent the color of the $\Phi_{i,j}$.
For the product of two symmetric representations S the situation is slightly
more complicated as one needs to represent each $\Phi_{u}$ where
$u=1,\cdots,\frac{1}{2}N(N+1)$ with two fundamental indices $i,j=1,\cdots,N$.
We can now write $A^{ij}_{kl}=v^{ij}w_{kl}$, which now has to be symmetric
under the transformations $i\leftrightarrow j$, $k\leftrightarrow l$.
Transforming between both representations can be done using the respective
Clebsch-Gordan coefficients Han:2009ya . The symmetricity representation
decomposes as

$\displaystyle A^{ij}_{kl}$
$\displaystyle=[\mathrm{\textbf{1}}]^{ij}_{kl}+[\mathrm{\textbf{A}}]^{ij}_{kl}+[\mathrm{\textbf{D}}]^{ij}_{kl}$
(86) $\displaystyle[\mathrm{\textbf{1}}]^{ij}_{kl}$
$\displaystyle=\frac{1}{N(N\\!+\\!1)}A^{mn}_{mn}\left(\delta^{i}_{k}\delta^{j}_{l}\\!+\\!\delta^{i}_{l}\delta^{j}_{k}\right)$
$\displaystyle[\mathrm{\textbf{A}}]^{ij}_{kl}$
$\displaystyle=\frac{1}{N\\!+\\!2}\left[\delta^{i}_{k}A^{mj}_{ml}\\!+\\!\delta^{j}_{k}A^{mi}_{ml}\\!+\\!\delta^{j}_{l}A^{mi}_{mk}\\!+\\!\delta^{i}_{l}A^{mj}_{mk}\right]\\!-\\!\frac{2}{N(N\\!+\\!2)}A^{mn}_{mn}\left(\delta^{i}_{k}\delta^{j}_{l}\\!+\\!\delta^{i}_{l}\delta^{j}_{k}\right)$
$\displaystyle[\mathrm{\textbf{D}}]^{ij}_{kl}$
$\displaystyle=A^{ij}_{kl}\\!-\\!\frac{1}{N\\!+\\!2}\left[\delta^{i}_{k}A^{mj}_{ml}\\!+\\!\delta^{j}_{k}A^{mi}_{ml}\\!+\\!\delta^{j}_{l}A^{mi}_{mk}\\!+\\!\delta^{i}_{l}A^{mj}_{mk}\right]\\!+\\!\frac{1}{(N\\!+\\!1)(N\\!+\\!2)}A^{mn}_{mn}\left(\delta^{i}_{k}\delta^{j}_{l}\\!+\\!\delta^{i}_{l}\delta^{j}_{k}\right).$

For the adjoint we can write $A^{ij}_{kl}=v^{i}_{k}w^{j}_{l}$ with
$A^{mj}_{ml}=A^{im}_{km}=0$. Again we can transform to adjoint indices
$a=1,\cdots,N^{2}-1$ by using the respective Clebsch-Gordan coefficients which
are obtained directly from the generators of $SU(N)$. The adjoint decomposes
as

$\displaystyle A^{ij}_{kl}=$
$\displaystyle[\mathrm{\textbf{1}}_{\textbf{S}}]^{ij}_{kl}+[\mathrm{\textbf{A}}_{\textbf{A}}]^{ij}_{kl}+[\mathrm{\textbf{A}}_{\textbf{S}}]^{ij}_{kl}+[\mathrm{\textbf{B}}_{\textbf{S}}]^{ij}_{kl}+[\mathrm{\textbf{C}}_{\textbf{A}}]^{ij}_{kl}+[\overline{\mathrm{\textbf{C}}}_{\textbf{A}}]^{ij}_{kl}+[\mathrm{\textbf{D}}_{\textbf{S}}]^{ij}_{kl}$
(87) $\displaystyle[\mathrm{\textbf{1}}_{\textbf{S}}]^{ij}_{kl}=$
$\displaystyle\frac{1}{N\\!-\\!N^{3}}A^{mn}_{nm}\left(\delta^{i}_{k}\delta^{j}_{l}-N\delta^{i}_{l}\delta^{j}_{k}\right)$
$\displaystyle[\mathrm{\textbf{A}}_{\textbf{A}}]^{ij}_{kl}=$
$\displaystyle\frac{1}{2N}\left[-\delta^{i}_{l}(A^{jm}_{mk}-A^{mj}_{km})+\delta^{j}_{k}(A^{im}_{ml}-A^{mi}_{lm})\right]$
$\displaystyle[\mathrm{\textbf{A}}_{\textbf{S}}]^{ij}_{kl}=$
$\displaystyle\frac{4}{N(N^{2}\\!-\\!4)}A^{mn}_{nm}\left(\delta^{i}_{k}\delta^{j}_{l}-\frac{N}{2}\delta^{i}_{l}\delta^{j}_{k}\right)$
$\displaystyle+\\!\frac{1}{4\\!-\\!N^{2}}\\!\left[\delta^{i}_{k}(A^{jm}_{ml}\\!+\\!A^{mj}_{lm})\\!-\\!\frac{N}{2}\delta^{i}_{l}(A^{jm}_{mk}\\!+\\!A^{mj}_{km})\\!-\\!\frac{N}{2}\delta^{j}_{k}(A^{im}_{ml}\\!+\\!A^{mi}_{lm})\\!+\\!\delta^{j}_{l}(A^{im}_{mk}\\!+\\!A^{mi}_{km})\right]$
$\displaystyle[\mathrm{\textbf{B}}_{\textbf{S}}]^{ij}_{kl}=$
$\displaystyle\frac{1}{4}(A^{ij}_{kl}-A^{ji}_{kl}-A^{ij}_{lk}+A^{ji}_{lk})-\frac{1}{2(N^{2}\\!-\\!3N\\!+\\!2)}A^{mn}_{nm}\left(\delta^{i}_{k}\delta^{j}_{l}-\delta^{i}_{l}\delta^{j}_{k}\right)$
$\displaystyle+\\!\frac{1}{4(N\\!-\\!2)}\\!\left[\delta^{i}_{k}(A^{jm}_{ml}\\!+\\!A^{mj}_{lm})\\!-\\!\delta^{i}_{l}(A^{jm}_{mk}\\!+\\!A^{mj}_{km})\\!-\\!\delta^{j}_{k}(A^{im}_{ml}\\!+\\!A^{mi}_{lm})\\!+\\!\delta^{j}_{l}(A^{im}_{mk}\\!+\\!A^{mi}_{km})\right]$
$\displaystyle[\mathrm{\textbf{C}}_{\textbf{A}}]^{ij}_{kl}=$
$\displaystyle\frac{1}{4}(A^{ij}_{kl}+A^{ji}_{kl}-A^{ij}_{lk}-A^{ji}_{lk})$
$\displaystyle-\\!\frac{1}{4N}\\!\left[\delta^{i}_{k}(A^{jm}_{ml}\\!-\\!A^{mj}_{lm})\\!-\\!\delta^{i}_{l}(A^{jm}_{mk}\\!-\\!A^{mj}_{km})\\!+\\!\delta^{j}_{k}(A^{im}_{ml}\\!-\\!A^{mi}_{lm})\\!-\\!\delta^{j}_{l}(A^{im}_{mk}\\!-\\!A^{mi}_{km})\right]$
$\displaystyle[\overline{\mathrm{\textbf{C}}}_{\textbf{A}}]^{ij}_{kl}=$
$\displaystyle\frac{1}{4}(A^{ij}_{kl}-A^{ji}_{kl}+A^{ij}_{lk}-A^{ji}_{lk})$
$\displaystyle+\\!\frac{1}{4N}\\!\left[\delta^{i}_{k}(A^{jm}_{ml}\\!-\\!A^{mj}_{lm})\\!+\\!\delta^{i}_{l}(A^{jm}_{mk}\\!-\\!A^{mj}_{km})\\!-\\!\delta^{j}_{k}(A^{im}_{ml}\\!-\\!A^{mi}_{lm})\\!-\\!\delta^{j}_{l}(A^{im}_{mk}\\!-\\!A^{mi}_{km})\right]$
$\displaystyle[\mathrm{\textbf{D}}_{\textbf{S}}]^{ij}_{kl}=$
$\displaystyle\frac{1}{4}(A^{ij}_{kl}+A^{ji}_{kl}+A^{ij}_{lk}+A^{ji}_{lk})+\frac{1}{2(N^{2}\\!+\\!3N\\!+\\!2)}A^{mn}_{nm}\left(\delta^{i}_{k}\delta^{j}_{l}+\delta^{i}_{l}\delta^{j}_{k}\right)$
$\displaystyle-\\!\frac{1}{4(N\\!+\\!2)}\\!\left[\delta^{i}_{k}(A^{jm}_{ml}\\!+\\!A^{mj}_{lm})\\!+\\!\delta^{i}_{l}(A^{jm}_{mk}\\!+\\!A^{mj}_{km})\\!+\\!\delta^{j}_{k}(A^{im}_{ml}\\!+\\!A^{mi}_{lm})\\!+\\!\delta^{j}_{l}(A^{im}_{mk}\\!+\\!A^{mi}_{km})\right].$

The above decomposition applies for $N\geq 4$, however, in the case of QCD
with $N=3$ the representation $\mathrm{\textbf{B}}_{\textbf{S}}$ does not
appear and the symmetric adjoint representation for $SU(3)$ is given by
$[\mathrm{\textbf{8}}_{\textbf{S}}]^{ij}_{kl}=[\mathrm{\textbf{A}}_{\textbf{S}}]^{ij}_{kl}+[\mathrm{\textbf{B}}_{\textbf{S}}]^{ij}_{kl}$.
This concludes the decomposition of the amplitudes considered in sections 3.2
and 5.1.

###  B.2 Decuplet annihilation

It is possible to imagine dark sectors with exotic and large representations
of $SU(3)$. Out of these representations the 10, 15 and 27 will have
annihilations directly into two Standard Model particles, either to two gluons
or a quark gluon pair. Although building these models is challenging, large
color representations are associated to higher annihilation cross sections
compared to the models we study — see equation (61). The large annihilation
rate leads to an enhanced depletion of the dark matter relic abundance or
equivalently allows for larger mass splittings between the dark matter and its
coannihilation partner. As an example we consider the 10 for which we have the
following color decomposition

$\mathrm{\textbf{10}}\otimes\overline{\mathrm{\textbf{10}}}=\mathrm{\textbf{1}}\oplus\mathrm{\textbf{8}}\oplus\mathrm{\textbf{27}}\oplus\mathrm{\textbf{64}}.$
(88)

By virtue of equation (40) while inserting the quadratic Casimir invariants
from equation (42) we decompose the QCD potential as

$V_{\mathrm{\textbf{10}}\otimes\overline{\mathrm{\textbf{10}}}}=\frac{\alpha_{s}}{r}\left\\{\begin{aligned}
-6&\quad(\mathrm{\textbf{1}})\\\ -\frac{9}{2}&\quad(\mathrm{\textbf{8}})\\\
-2&\quad(\mathrm{\textbf{27}})\\\
+\frac{3}{2}&\quad(\mathrm{\textbf{64}})\end{aligned}\right..$ (89)

To decompose the total cross section we write the color part of the amplitude
in tensor notation as $A^{ijk}_{lmn}=v^{ijk}w_{lmn}$ with full symmetricity in
the upper and lower components separately. After doing the calculation we find

$\displaystyle{[\mathrm{\textbf{1}}]}^{ijk}_{lmn}$
$\displaystyle=\frac{1}{60}D^{ijk}_{lmn}$ (90)
$\displaystyle{[\mathrm{\textbf{8}}]}^{ijk}_{lmn}$
$\displaystyle=\frac{1}{30}C^{ijk}_{lmn}-\frac{1}{30}D^{ijk}_{lmn}$
$\displaystyle{[\mathrm{\textbf{27}}]}^{ijk}_{lmn}$
$\displaystyle=\frac{1}{7}B^{ijk}_{lmn}-\frac{2}{35}C^{ijk}_{lmn}+\frac{3}{140}D^{ijk}_{lmn}$
$\displaystyle{[\mathrm{\textbf{64}}]}^{ijk}_{lmn}$
$\displaystyle=A^{ijk}_{lmn}-\frac{1}{7}B^{ijk}_{lmn}+\frac{1}{42}C^{ijk}_{lmn}-\frac{1}{210}D^{ijk}_{lmn}.$

In these equations we used

$\displaystyle B^{ijk}_{lmn}$
$\displaystyle=\delta^{i}_{l}A^{pjk}_{pmn}+\cdots$ (91) $\displaystyle
C^{ijk}_{lmn}$
$\displaystyle=\delta^{i}_{l}\delta^{j}_{m}A^{pqk}_{pqn}+\cdots$
$\displaystyle D^{ijk}_{lmn}$
$\displaystyle=\delta^{i}_{l}\delta^{j}_{m}\delta^{k}_{n}A^{pqr}_{pqr}+\cdots,$

where the dots represent all symmetric combinations in the upper and lower
indices. $B^{ijk}_{lmn}$ has nine terms, $C^{ijk}_{lmn}$ has eighteen terms
and $D^{ijk}_{lmn}$ has six terms. By inserting equation (47) into equation
(90), we then obtain the following decomposition for the
$\mathrm{\textbf{10}}\otimes\overline{\mathrm{\textbf{10}}}\to g^{a}\,g^{b}$
process

$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{10}}\otimes\overline{\mathrm{\textbf{10}}}}\big{|}^{2}$
$\displaystyle=7\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}]\big{|}^{2}+\frac{35}{9}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2}+\frac{5}{3}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{27}}]\big{|}^{2}$
$\displaystyle\quad\mathrm{even}\,\,l+s$ (92)
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{10}}\otimes\overline{\mathrm{\textbf{10}}}}\big{|}^{2}$
$\displaystyle=\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2}$
$\displaystyle\quad\mathrm{odd}\,\,l+s$ .

As expected the amplitudes for odd $l+s$ only involve the color-octet channel.
Moreover, for even $l+s$ we observe no decomposition into the 64 as this
representation does not appear in the color product of two gluons. From this
result and equation (89) the Sommerfeld corrections are obtained equivalently
to equation (56) as

$\sigma^{(S)}_{\mathrm{\textbf{10}}\otimes\overline{\mathrm{\textbf{10}}}\to
g\,g}=\begin{cases}\frac{1}{7}\sigma^{(S)}_{C}\left[6\alpha_{s}\right]+\frac{9}{35}\sigma^{(S)}_{C}\left[\frac{9\alpha_{s}}{2}\right]+\frac{3}{5}\sigma^{(S)}_{C}\left[2\alpha_{s}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[\frac{9\alpha_{s}}{2}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}.$
(93)

We observe here that in contrast to the correction factors for the triplet,
sextet and octet in equation (56), the decuplet has positive coupling
strengths for all of the Coulomb potentials. This implies an even larger
enhancement of the annihilation cross sections and strengthens the effect of
the larger Casimir values.

###  B.3 Triplet–octet annihilation

An interesting possible scenario is a dark sector with two colored particles
close in mass to the dark matter particle. For example, one could consider a
model with a triplet and an octet of $SU(3)$. In addition to the self-
annihilation of each of these particles through the strong interaction, the 3
and the 8 could coannihilate to a quark and a gluon. The color of the initial
state can be decomposed as

$\mathrm{\textbf{3}}\otimes\mathrm{\textbf{8}}=\mathrm{\textbf{3}}\oplus\overline{\mathrm{\textbf{6}}}\oplus\mathrm{\textbf{15}}.$
(94)

We then use equation (40) and (42) to decompose the QCD potential which gives

$V_{\mathrm{\textbf{3}}\otimes\mathrm{\textbf{8}}}=\frac{\alpha_{s}}{r}\left\\{\begin{aligned}
-\frac{3}{2}&\quad(\mathrm{\textbf{3}})\\\
-\frac{1}{2}&\quad(\overline{\mathrm{\textbf{6}}})\\\
+\frac{1}{2}&\quad(\mathrm{\textbf{15}})\end{aligned}\right..$ (95)

To decompose the total cross section we write the color part of the amplitude
in tensor notation as $A^{ij}_{k}=v^{i}w^{j}_{k}$ with the condition
$A^{im}_{m}=0$. We then find

$\displaystyle{[\mathrm{\textbf{3}}]}^{ij}_{k}$
$\displaystyle=\frac{3}{8}\delta^{i}_{k}A^{mj}_{m}-\frac{1}{8}\delta^{j}_{k}A^{mi}_{m}$
(96) $\displaystyle{[\overline{\mathrm{\textbf{6}}}]}^{ij}_{k}$
$\displaystyle=\frac{1}{2}(A^{ij}_{k}-A^{ji}_{k})+\frac{1}{4}(\delta^{j}_{k}A^{mi}_{m}-\delta^{i}_{k}A^{mj}_{m})$
$\displaystyle{[\mathrm{\textbf{15}}]}^{ij}_{k}$
$\displaystyle=\frac{1}{2}(A^{ij}_{k}+A^{ji}_{k})-\frac{1}{8}(\delta^{j}_{k}A^{mi}_{m}+\delta^{i}_{k}A^{mj}_{m}).$

In this model we assume a coupling between new particles that transform under
the 3 and the 8 and a Standard Model quark — such as the squark-gluino
coupling in supersymmetry. The
$\mathrm{\textbf{3}}\otimes\mathrm{\textbf{8}}\to q\,g$ coannihilation process
can then occur through either an $s$-channel quark, a $t$-channel 3 or a
$t$-channel 8. These diagrams have different color structures and their
relative strength determines the decomposition over the three different color
channels. In contrast to the processes considered in the main part an $l+s$
symmetry is not applicable to this process since neither the initial nor the
final state involve pairs of identical particles. Therefore a decomposition
like

$\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{3}}\otimes\mathrm{\textbf{8}}}\big{|}^{2}=\alpha\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{3}}]\big{|}^{2}=\beta\sum_{\mathrm{color}}\big{|}[\overline{\mathrm{\textbf{6}}}]\big{|}^{2}=\gamma\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{15}}]\big{|}^{2}.$
(97)

will have momentum-dependent factors $\alpha$, $\beta$ and $\gamma$, rendering
the calculation of Sommerfeld corrections to be more cumbersome.

A possible approach is based on the fact that each of the three diagrams
contributing to the amplitude has a fixed color structure. When squaring the
amplitude, the squares of the contributions of each diagram as well the
interference terms will also have a well-defined color structure when
considered separately. Applying the recipe described in section 3.3 to each of
these terms will give the total analytic Sommerfeld-corrected cross section at
a given order in the partial wave expansion. For convenience we present here
the color decompositions of each combination of channels:

$\displaystyle\alpha=1,\qquad$ $\displaystyle\beta=0,\qquad$
$\displaystyle\gamma=0$ $\displaystyle\qquad s\textrm{-channel squared}$ (98)
$\displaystyle\alpha=64,\qquad$ $\displaystyle\beta=\frac{32}{9},\qquad$
$\displaystyle\gamma=\frac{64}{45}$ $\displaystyle\qquad t_{3}\textrm{-channel
squared}$ $\displaystyle\alpha=\frac{16}{9},\qquad$
$\displaystyle\beta=8,\qquad$ $\displaystyle\gamma=\frac{16}{5}$
$\displaystyle\qquad t_{8}\textrm{-channel squared}$
$\displaystyle\alpha=1,\qquad$ $\displaystyle\beta=0,\qquad$
$\displaystyle\gamma=0$ $\displaystyle\qquad s\textrm{-channel interference}$
$\displaystyle\alpha=8,\qquad$ $\displaystyle\beta=4,\qquad$
$\displaystyle\gamma=\frac{8}{5}$ $\displaystyle\qquad t_{3/8}\textrm{-channel
interference}.$

###  B.4 Triplet–triplet annihilation

Here, we consider a model involving a scalar dark matter particle that couples
to a Standard Model quark and a new vector-like quark ($\psi$) which is a
triplet under color Giacchino:2015hvk . In this case, the annihilation of
$\psi$ into quark pairs can occur through two independent processes, namely
$\psi\,\bar{\psi}\to q\,\bar{q}$ and $\psi\,\psi\to q\,q$ plus its conjugate.
Annihilation into a quark anti-quark pair can occur through either an
$s$-channel gluon or a $t$-channel scalar (the dark matter). $\psi$
annihilation to identical quarks occurs either in the $t$ and $u$-channel
through dark matter exchange.

First we discuss the effect of the $t$-channel dark matter exchange to
$\psi\,\bar{\psi}\to q\,\bar{q}$. The corresponding potential has been derived
in (43) and the resulting Sommerfeld correction for the $s$-channel gluon
exchange has been derived in section 3.2 and presented in equation (56).
Including the new $t$-channel diagrams leads to a more complex color
structure. As for the triplet-octet model, since neither the initial state nor
the final state involve identical particles, the $(l,s)$ components of the
amplitude cannot be constrained by symmetry arguments. We have to adopt the
same strategy as in section B.3 and observe that the $t$-channel amplitude
decomposes as

$\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{3}}\otimes\bar{\mathrm{\textbf{3}}}}^{t\textrm{-channel}}\big{|}^{2}=9\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{1}}]\big{|}^{2}=\frac{9}{8}\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{8}}]\big{|}^{2},$
(99)

whereas interference between the $s$ and $t$-channel only occurs when the
initial state is in the octet representation.

The situation changes for the second process $\psi\,\psi\to q\,q$. Since the
quarks in the final state are identical, contributions from each $(l,s)$ state
are constrained by symmetry. The color part of this amplitude can be
decomposed as

$\mathrm{\textbf{3}}\otimes\mathrm{\textbf{3}}=\overline{\mathrm{\textbf{3}}}\oplus\mathrm{\textbf{6}}.$
(100)

We then use equations (40) and (42) to decompose the QCD potential which gives

$V_{\mathrm{\textbf{3}}\otimes\mathrm{\textbf{3}}}=\frac{\alpha_{s}}{r}\left\\{\begin{aligned}
-\frac{2}{3}&\quad(\mathrm{\textbf{3}})\\\
+\frac{1}{3}&\quad(\overline{\mathrm{\textbf{6}}})\end{aligned}\right..$ (101)

To decompose the total cross section we write the color part of the initial
states in tensor notation as $A^{ij}=v^{i}w^{j}$ and find

$\displaystyle{[\overline{\mathrm{\textbf{3}}}]}^{ij}$
$\displaystyle=\frac{1}{2}(A^{ij}-A^{ji})$ (102)
$\displaystyle{[\mathrm{\textbf{6}}]}^{ij}$
$\displaystyle=\frac{1}{2}(A^{ij}+A^{ji}).$

The total color structure for this process can be written as the sum of the
$t$-channel and $u$-channel contributions, namely
$A^{ij}_{kl}=\alpha\delta^{i}_{k}\delta^{j}_{l}+\beta\delta^{j}_{k}\delta^{i}_{l}$.
The CP symmetry condition akin to that of equation (48) imposes that
$\alpha=\beta$ for $l+s$ even and $\alpha=-\beta$ for $l+s$ odd. Inserting
this information into equation (102) gives the following color decomposition

$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{3}}\otimes\mathrm{\textbf{3}}}\big{|}^{2}$
$\displaystyle=\sum_{\mathrm{color}}\big{|}[\mathrm{\textbf{6}}]\big{|}^{2}$
$\displaystyle\quad\mathrm{even}\,\,l+s$ (103)
$\displaystyle\sum_{\mathrm{color}}\big{|}A_{\mathrm{\textbf{3}}\otimes\mathrm{\textbf{3}}}\big{|}^{2}$
$\displaystyle=\sum_{\mathrm{color}}\big{|}[\overline{\mathrm{\textbf{3}}}]\big{|}^{2}$
$\displaystyle\quad\mathrm{odd}\,\,l+s$ ,

which, using equation (101), leads to

$\sigma^{(S)}_{\mathrm{\textbf{3}}\otimes\mathrm{\textbf{3}}\to
q\,q}=\begin{cases}\sigma^{(S)}_{C}\left[-\frac{\alpha_{s}}{3}\right]&\mathrm{even}\,\,l+s\\\
\sigma^{(S)}_{C}\left[\frac{2\alpha_{s}}{3}\right]&\phantom{!}\mathrm{odd}\,\,l+s\end{cases}$
(104)

which, beyond the $s$-wave, contrasts with the result derived in
Giacchino:2015hvk  for a similar simplified model.

## References

  * (1) A. Sommerfeld, Über die Beugung und Bremsung der Elektronen, Ann. Phys. 403 (1931) 257–330.
  * (2) J. Hisano, S. Matsumoto, M. M. Nojiri, and O. Saito, Non-perturbative effect on dark matter annihilation and gamma ray signature from galactic center, Phys. Rev. D71 (2005) 063528, [hep-ph/0412403].
  * (3) N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, A Theory of Dark Matter, Phys. Rev. D79 (2009) 015014, [arXiv:0810.0713].
  * (4) M. Lattanzi and J. I. Silk, Can the WIMP annihilation boost factor be boosted by the Sommerfeld enhancement?, Phys. Rev. D79 (2009) 083523, [arXiv:0812.0360].
  * (5) T. Cohen, M. Lisanti, A. Pierce, and T. R. Slatyer, Wino Dark Matter Under Siege, JCAP 1310 (2013) 061, [arXiv:1307.4082].
  * (6) J. Fan and M. Reece, In Wino Veritas? Indirect Searches Shed Light on Neutralino Dark Matter, JHEP 10 (2013) 124, [arXiv:1307.4400].
  * (7) A. Hryczuk, I. Cholis, R. Iengo, M. Tavakoli, and P. Ullio, Indirect Detection Analysis: Wino Dark Matter Case Study, JCAP 1407 (2014) 031, [arXiv:1401.6212].
  * (8) M. Beneke, A. Bharucha, A. Hryczuk, S. Recksiegel, and P. Ruiz-Femenia, The last refuge of mixed wino-Higgsino dark matter, JHEP 01 (2017) 002, [arXiv:1611.00804].
  * (9) J. Hisano, S. Matsumoto, M. Nagai, O. Saito, and M. Senami, Non-perturbative effect on thermal relic abundance of dark matter, Phys. Lett. B646 (2007) 34–38, [hep-ph/0610249].
  * (10) A. Hryczuk, R. Iengo, and P. Ullio, Relic densities including Sommerfeld enhancements in the MSSM, JHEP 03 (2011) 069, [arXiv:1010.2172].
  * (11) M. Beneke, C. Hellmann, and P. Ruiz-Femenia, Non-relativistic pair annihilation of nearly mass degenerate neutralinos and charginos III. Computation of the Sommerfeld enhancements, JHEP 05 (2015) 115, [arXiv:1411.6924].
  * (12) M. Beneke, A. Bharucha, F. Dighera, C. Hellmann, A. Hryczuk, S. Recksiegel, and P. Ruiz-Femenia, Relic density of wino-like dark matter in the MSSM, JHEP 03 (2016) 119, [arXiv:1601.04718].
  * (13) Planck Collaboration, P. A. R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, arXiv:1502.01589.
  * (14) M. J. Baker et al., The Coannihilation Codex, JHEP 12 (2015) 120, [arXiv:1510.03434].
  * (15) M. Buschmann, S. El Hedri, A. Kaminska, J. Liu, M. de Vries, X.-P. Wang, F. Yu, and J. Zurita, Hunting for dark matter coannihilation by mixing dijet resonances and missing transverse energy, JHEP 09 (2016) 033, [arXiv:1605.08056].
  * (16) A. Strumia, Sommerfeld corrections to type-II and III leptogenesis, Nucl. Phys. B809 (2009) 308–317, [arXiv:0806.1630].
  * (17) A. De Simone, G. F. Giudice, and A. Strumia, Benchmarks for Dark Matter Searches at the LHC, JHEP 06 (2014) 081, [arXiv:1402.6287].
  * (18) J. Ellis, F. Luo, and K. A. Olive, Gluino Coannihilation Revisited, JHEP 09 (2015) 127, [arXiv:1503.07142].
  * (19) A. Ibarra, A. Pierce, N. R. Shah, and S. Vogl, Anatomy of Coannihilation with a Scalar Top Partner, Phys. Rev. D91 (2015), no. 9 095018, [arXiv:1501.03164].
  * (20) J. Harz, B. Herrmann, M. Klasen, K. Kovařík, and M. Meinecke, SUSY-QCD corrections to stop annihilation into electroweak final states including Coulomb enhancement effects, Phys. Rev. D91 (2015), no. 3 034012, [arXiv:1410.8063].
  * (21) S. Cassel, Sommerfeld factor for arbitrary partial wave processes, J. Phys. G37 (2010) 105009, [arXiv:0903.5307].
  * (22) S. Cassel, D. M. Ghilencea, and G. G. Ross, Electroweak and Dark Matter Constraints on a Z-prime in Models with a Hidden Valley, Nucl. Phys. B827 (2010) 256–280, [arXiv:0903.1118].
  * (23) R. Iengo, Sommerfeld enhancement: General results from field theory diagrams, JHEP 05 (2009) 024, [arXiv:0902.0688].
  * (24) T. Gherghetta, M. Nguyen, and Z. Thomas, Neutral Naturalness with Bifundamental Gluinos, Phys. Rev. D94 (2016), no. 11 115008, [arXiv:1610.00342].
  * (25) M. A. Buen-Abad, G. Marques-Tavares, and M. Schmaltz, Non-Abelian dark matter and dark radiation, Phys. Rev. D92 (2015), no. 2 023531, [arXiv:1505.03542].
  * (26) P. Ko and Y. Tang, Residual Non-Abelian Dark Matter and Dark Radiation, Phys. Lett. B768 (2017) 12–17, [arXiv:1609.02307].
  * (27) J. L. Feng and Y. Shadmi, WIMPless Dark Matter from Non-Abelian Hidden Sectors with Anomaly-Mediated Supersymmetry Breaking, Phys. Rev. D83 (2011) 095011, [arXiv:1102.0282].
  * (28) K. K. Boddy, J. L. Feng, M. Kaplinghat, Y. Shadmi, and T. M. P. Tait, Strongly interacting dark matter: Self-interactions and keV lines, Phys. Rev. D90 (2014), no. 9 095016, [arXiv:1408.6532].
  * (29) K. K. Boddy, J. L. Feng, M. Kaplinghat, and T. M. P. Tait, Self-Interacting Dark Matter from a Non-Abelian Hidden Sector, Phys. Rev. D89 (2014), no. 11 115017, [arXiv:1402.3629].
  * (30) K. Harigaya, M. Ibe, K. Kaneta, W. Nakano, and M. Suzuki, Thermal Relic Dark Matter Beyond the Unitarity Limit, JHEP 08 (2016) 151, [arXiv:1606.00159].
  * (31) A. Soni, H. Xiao, and Y. Zhang, A Cosmic Selection Rule for Glueball Dark Matter Relic Density, arXiv:1704.02347.
  * (32) S. El Hedri, A. Kaminska, M. de Vries, and J. Zurita, Simplified Phenomenology for Colored Dark Sectors, JHEP 04 (2017) 118, [arXiv:1703.00452].
  * (33) L. Visinelli and P. Gondolo, An integral equation for distorted wave amplitudes, arXiv:1007.2903.
  * (34) W. Fischler, Quark - anti-Quark Potential in QCD, Nucl. Phys. B129 (1977) 157–174.
  * (35) Y. Schröder, The Static potential in QCD to two loops, Phys. Lett. B447 (1999) 321–326, [hep-ph/9812205].
  * (36) M. R. Kauth, J. H. Kuhn, P. Marquard, and M. Steinhauser, Gluino Pair Production at the LHC: The Threshold, Nucl. Phys. B857 (2012) 28–64, [arXiv:1108.0361].
  * (37) F. Giacchino, A. Ibarra, L. Lopez Honorez, M. H. G. Tytgat, and S. Wild, Signatures from Scalar Dark Matter with a Vector-like Quark Mediator, JCAP 1602 (2016), no. 02 002, [arXiv:1511.04452].
  * (38) H. Georgi, Lie algebras in particle physics, Front. Phys. 54 (1999) 1–320.
  * (39) S. P. Liew and F. Luo, Effects of QCD bound states on dark matter relic abundance, JHEP 02 (2017) 091, [arXiv:1611.08133].
  * (40) C. F. Berger, L. Covi, S. Kraml, and F. Palorini, The Number density of a charged relic, JCAP 0810 (2008) 005, [arXiv:0807.0211].
  * (41) H. Baer, K.-m. Cheung, and J. F. Gunion, A Heavy gluino as the lightest supersymmetric particle, Phys. Rev. D59 (1999) 075002, [hep-ph/9806361].
  * (42) J. Blumlein, E. Boos, and A. Kryukov, Leptoquark pair production in hadronic interactions, Z. Phys. C76 (1997) 137–153, [hep-ph/9610408].
  * (43) J. L. Hewett, T. G. Rizzo, S. Pakvasa, H. E. Haber, and A. Pomarol, Vector leptoquark production at hadron colliders, in Workshop on Physics at Current Accelerators and the Supercollider Argonne, Illinois, June 2-5, 1993, 1993.  hep-ph/9310361.
  * (44) S. El Hedri, A. Kaminska, and M. de Vries, “Mathematica Notebook for Analytic Sommerfeld Corrections.”  GitHub repository: https://github.com/MDT-Maikel/Sommerfeld.
  * (45) G. Bélanger, F. Boudjema, A. Pukhov, and A. Semenov, micrOMEGAs4.1: two dark matter candidates, Comput. Phys. Commun. 192 (2015) 322–329, [arXiv:1407.6129].
  * (46) D. Barducci, G. Belanger, J. Bernon, F. Boudjema, J. Da Silva, S. Kraml, U. Laa, and A. Pukhov, Collider limits on new physics within micrOMEGAs, arXiv:1606.03834.
  * (47) Particle Data Group Collaboration, C. Patrignani et al., Review of Particle Physics, Chin. Phys. C40 (2016), no. 10 100001. Chapter 9. Quantum Chromodynamics.
  * (48) G. M. Prosperi, M. Raciti, and C. Simolo, On the running coupling constant in QCD, Prog. Part. Nucl. Phys. 58 (2007) 387–438, [hep-ph/0607209].
  * (49) J. E. Juknevich, D. Melnikov, and M. J. Strassler, A Pure-Glue Hidden Valley I. States and Decays, JHEP 07 (2009) 055, [arXiv:0903.0883].
  * (50) J. E. Juknevich, Pure-glue hidden valleys through the Higgs portal, JHEP 08 (2010) 121, [arXiv:0911.5616].
  * (51) S. Jeon and R. Venugopalan, Random walks of partons in SU(N(c)) and classical representations of color charges in QCD at small x, Phys. Rev. D70 (2004) 105012, [hep-ph/0406169].
  * (52) S. Okubo, Casimir Invariants and Vector Operators in Simple Lie Algebra, J. Math. Phys. 18 (1977) 2382.
  * (53) Z. Zhang, Multi-Sommerfeld enhancement in dark sector, Phys. Lett. B734 (2014) 188–192, [arXiv:1307.2206].
  * (54) O. Just, The Partial Wave Formalism and its Application to Neutralino Dark Matter, Master’s thesis, 2008.  http://wwwmpa.mpa-garching.mpg.de/~ojust/ojust_files/Just_diplomarbeit.pdf.
  * (55) S. U. Chung, Spin Formalisms, CERN-71-08 (1971). https://suchung.web.cern.ch/suchung/spinfm1.pdf.
  * (56) R. Feger and T. W. Kephart, LieART—A Mathematica application for Lie algebras and representation theory, Comput. Phys. Commun. 192 (2015) 166–195, [arXiv:1206.6379].
  * (57) M. Sjödahl, ColorMath - A package for color summed calculations in SU(Nc), Eur. Phys. J. C73 (2013), no. 2 2310, [arXiv:1211.2099].
  * (58) T. Han, I. Lewis, and T. McElmurry, QCD Corrections to Scalar Diquark Production at Hadron Colliders, JHEP 01 (2010) 123, [arXiv:0909.2666].