1 Introduction

Sommerfeld corrections [1] 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 [2,3,4]. In particular, the predicted enhancement of the annihilation rate in the galactic center for pure wino dark matter has allowed one to strongly restrict this supersymmetric scenario [5,6,7,8]. For multi-TeV weakly interacting dark matter, Sommerfeld corrections also typically lead to order one modifications of its relic density [9,10,11,12], often significantly weakening the upper bound on the dark matter mass derived from the Planck measurement [13].

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 [14, 15]—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 [16, 17]. 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 [16,17,18,19,20]. Notably, reference [17] 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 [21,22,23]. 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 the face of growing interest in the cosmological role of new gauge groups [24,25,26,27,28,29,30,31], 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 [32], 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 Sect. 2 we discuss the analytic derivation of Sommerfeld corrections to annihilation processes for arbitrary partial waves and with any momentum dependence. In Sect. 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 Sects. 2 and 3 are orthogonal and can be combined into a general prescription for the annihilation of colored particles. In Sect. 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 Sect. 5. We conclude in Sect. 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

$$\begin{aligned} {\mathcal {M}}(p, \theta , \phi ) = \sum _{l,m} F_{lm}(p) Y_{lm} (\theta ,\phi ), \end{aligned}$$
(2.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\)

$$\begin{aligned} F_{lm}(p) = \sum _{n\ge 0} \alpha _{lmn} \, p^{l+2n}. \end{aligned}$$
(2.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,

$$\begin{aligned} \sigma _{\mathrm {Sommerfeld}} = S \sigma _{\mathrm {perturbative}}. \end{aligned}$$
(2.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 [1]

$$\begin{aligned} S(x) = \frac{\pi x}{1 - \mathrm{e}^{-\pi x}}, \quad x = \frac{A}{\beta }. \end{aligned}$$
(2.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 formulas for the Sommerfeld-correction factors for higher waves have been computed in [21, 23] assuming the amplitude is proportional to \(p^l\) for the lth partial wave. This has been extended upon slightly in [22] allowing for a single term with a momentum dependence of \(p^{l+2n}\) with \(n \ge 0\). Here, we extend these results to a full expansion of the annihilation amplitude into orbital angular momentum and spin states (ls) 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 Eq. (2.2) into account. In the rest of this section we consider a Coulomb potential and do not make assumptions as regards the spin of the initial-state particles.

2.1 Partial wave expansion

For a given field \(\Phi \), the \(\Phi \, \overline{\Phi } \rightarrow \text {SM} \, \text {SM}\) amplitude can be expanded into orbital angular momentum and spin states (ls). 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 Sect. 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 [22], 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 one to obtain analytical formulas for the Sommerfeld-corrected matrix element \({\mathcal {M}}^{(S)}_{ls}\). Note that, since the different (ls) states are orthogonal to each other, the final cross section will be of the form

$$\begin{aligned} \sigma ^{(S)} \propto \sum _{l,s} \left| {\mathcal {M}}^{(S)}_{ls} \right| ^2. \end{aligned}$$
(2.5)

Another notable advantage of using a (ls) decomposition is that for identical particles the overall form of a given (ls) 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 } \rightarrow 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_{\mathbf{R }}, T^b_{\mathbf{R }} \right\} \) while the ones for states with odd \(l+s\) will be proportional to \(\left[ T^a_{\mathbf{R }}, T^b_{\mathbf{R }} \right] \) where \(T^a_{\mathbf{R }}\) is the generator for the representation \(\mathbf{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 Sect. 3. The same arguments apply to the case of SU(N), which will be discussed in Sect. 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:

$$\begin{aligned} {\mathcal {T}}_{fi} \left( p, \theta , \phi \right) {=} \langle p_f; 0 0; m_3 m_4 | T | p; \theta \phi ; m_1 m_2 \rangle . \end{aligned}$$
(2.6)

Here, p and \(p_f\) being the magnitudes of the momenta in the initial and final states, respectively. The information as regards the total spins \(s_{1,2,3,4}\) in the initial and final state is omitted here for compactness of notation. Further details as regards 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 asFootnote 1

$$\begin{aligned}&{\mathcal {T}}_{fi} \left( p, \theta , \phi \right) \nonumber \\&\quad =\sum _{l, l_z} \langle p_f; 0 0; m_3 m_4 | T | p; l l_z; m_1 m_2 \rangle \, Y^{l_z}_l(\theta , \phi ). \end{aligned}$$
(2.7)

A given \(| p_i; l l_z; m_1 m_2 \rangle \) state can be decomposed into \(| p; l l_z; s s_z \rangle \) states using Clebsch–Gordan coefficients,

$$\begin{aligned} | p; l l_z; m_1 m_2 \rangle = \sum _{s, s_z} \langle s_1 m_1 s_2 m_2 | s s_z \rangle | p; l l_z; s s_z \rangle , \end{aligned}$$
(2.8)

where \(s_1, s_2\) are the total spins of the incoming particles. The total amplitude can then be written as

$$\begin{aligned} {\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 | s s_z \rangle \nonumber \\&\times \, {\mathcal {M}}(p; l l_z; s s_z; m_3 m_4). \end{aligned}$$
(2.9)

The matrix element \({\mathcal {M}}(p; l l_z; s s_z; m_3 m_4)\) corresponds to the contribution of a single initial state \(| p_i; l l_z; s s_z \rangle \) to the total amplitude. Knowing \({\mathcal {T}}_{fi}\), this matrix element can be computed using

$$\begin{aligned} {\mathcal {M}}(p; l l_z; s s_z; m_3 m_4)\equiv & {} \langle p_f; 0 0; m_3 m_4 | T | p; l l_z; s m_s \rangle \nonumber \\= & {} \sum _{m_1, m_2} \langle s_1 m_1 s_2 m_2 | s s_z \rangle \nonumber \\&\times \int {\mathrm {d}} \Omega \, Y^{l_z \, *}_l(\theta , \phi ) {\mathcal {T}}_{fi} \left( p, \theta , \phi \right) .\nonumber \\ \end{aligned}$$
(2.10)

Since the (ls) components of the amplitude are orthogonal the total cross section is of the form

$$\begin{aligned} \sigma= & {} \frac{1}{64 \pi ^2 s} \frac{1}{\sqrt{1 - \frac{4m_\Phi ^2}{s}}} \dfrac{1}{d_\Phi ^2 d_{\mathbf{R }}^2}\nonumber \\&\times \sum _{m_3, m_4} \sum _{l, l_z} \sum _{s, s_z} \left| {\mathcal {M}}(p; l l_z; s s_z; m_3 m_4) \right| ^2, \end{aligned}$$
(2.11)

where \(d_\Phi \) is the number of degrees of freedom of the field \(\Phi \) and \(d_{\mathbf{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

$$\begin{aligned} {\mathcal {M}}(p; l l_z; s s_z; m_3 m_4) = \sum _{n\ge 0} \alpha ^{(m_3,m_4)}_{ll_zss_z,n} \, p^{l+2n}. \end{aligned}$$
(2.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 Eq. (2.12). This will be derived in the next section and the total Sommerfeld-corrected cross sections can then be obtained by the use of Eq. (2.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 [22]. 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 Fig. 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 [33]:

$$\begin{aligned} {\mathcal {M}}^{(S)}_{\beta \alpha } = {\mathcal {M}}^0_{\beta \alpha } + \int \mathrm{d}\gamma \frac{{\mathcal {M}}^{(S)}_{\beta \gamma } V_{\gamma \alpha }}{E_\alpha - E_\gamma + i \epsilon }, \end{aligned}$$
(2.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

$$\begin{aligned} 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 , \end{aligned}$$
(2.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

$$\begin{aligned} V_{\gamma \alpha }= & {} \langle m_1 m_2 | m_a m_b \rangle \langle k; \theta _k \phi _k | \hat{V} | p; \theta _p \phi _p \rangle \nonumber \\= & {} \delta _{m_1 m_a} \delta _{m_2 m_b} V(|{\mathbf {p}} - {\mathbf {k}}|). \end{aligned}$$
(2.15)
Fig. 1
figure 1

Sommerfeld ladder diagram for the annihilation of \(\Phi \) into Standard Model particles

For initial- and final-state spins \(\mathbf{m }_i = \left\{ m_1, m_2 \right\} \) and \(\mathbf{m }_f = \left\{ m_3, m_4 \right\} \), the Sommerfeld-corrected amplitude can be expressed as

$$\begin{aligned} {\mathcal {M}}^{(S)}_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {p}})= & {} {\mathcal {M}}^0_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {p}})\nonumber \\&+ \int \frac{{\mathrm {d}}^3k}{(2 \pi )^3} \frac{{\mathcal {M}}^{(S)}_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {k}}) V(|{\mathbf {p}} - {\mathbf {k}}|)}{E_\alpha - E_\gamma + i \epsilon }. \end{aligned}$$
(2.16)

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

$$\begin{aligned} (E_\alpha - E_\gamma )^{-1} \approx {\mathcal {E}} - \frac{{\mathbf {k}}^2}{2 \mu }, \end{aligned}$$
(2.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 } \ne {\mathcal {E}}\), we can define [23]

$$\begin{aligned} \Phi _{\mathbf{m }_f \mathbf{m }_i}(\widetilde{{\mathbf {p}}}) = \frac{{\mathcal {M}}^{(S)}_{\mathbf{m }_f \mathbf{m }_i}(\widetilde{{\mathbf {p}}})}{\frac{\widetilde{{\mathbf {p}}}^2}{2 \mu } - {\mathcal {E}}}, \end{aligned}$$
(2.18)

which verifies

$$\begin{aligned}&\left( \frac{\widetilde{{\mathbf {p}}}^2}{2 \mu } -{\mathcal {E}}\right) \Phi _{\mathbf{m }_f \mathbf{m }_i} (\widetilde{{\mathbf {p}}})\nonumber \\&\quad ={\mathcal {M}}^0_{\mathbf{m }_f \mathbf{m }_i} (\widetilde{{\mathbf {p}}}) -\int \frac{{\mathrm {d}}^3 k}{(2 \pi )^3} \Phi _{\mathbf{m }_f \mathbf{m }_i}({\mathbf {k}}) V(|\widetilde{{\mathbf {p}}} - {\mathbf {k}}|). \end{aligned}$$
(2.19)

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

$$\begin{aligned} \left( \frac{-\mathbf {\nabla }^2}{2 \mu } + V(r) - {\mathcal {E}} \right) \widetilde{\Phi }_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {r}}) = U^0_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {r}}), \end{aligned}$$
(2.20)

where

$$\begin{aligned} U^0_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {r}}) = \int \frac{\mathrm{d}^3 q}{(2 \pi )^3} \mathrm{e}^{i {\mathbf {r}} \cdot {\mathbf {q}}} {\mathcal {M}}^0_{\mathbf{m }_f \mathbf{m }_i} ({\mathbf {q}}). \end{aligned}$$
(2.21)

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

$$\begin{aligned}&{\mathcal {M}}^{(S)}_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {p}}) = \lim _{\widetilde{{\mathbf {p}}}\rightarrow {\mathbf {p}}} \left( \frac{\widetilde{{\mathbf {p}}}^2}{2 \mu }- {\mathcal {E}} \right) \Phi _{\mathbf{m }_f \mathbf{m }_i}(\widetilde{{\mathbf {p}}}) \nonumber \\&\quad {\mathrm {with}}\ \frac{{\mathbf {p}}^2}{2 \mu } = {\mathcal {E}}, \end{aligned}$$
(2.22)

which leads to [21, 23]

$$\begin{aligned} {\mathcal {M}}^{(S)}_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {p}}) = \int \frac{{\mathrm {d}}^3 q}{(2 \pi )^3} \, {\mathcal {M}}^0_{\mathbf{m }_f \mathbf{m }_i}({\mathbf {q}}) \, \phi _{{\mathbf {p}}}({\mathbf {q}}), \end{aligned}$$
(2.23)

where \(\phi _{{\mathbf {p}}}({\mathbf {q}})\) obeys the traditional Schrödinger equation in position space,

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

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

$$\begin{aligned} \phi _{{\mathbf {p}}}({\mathbf {q}}) = \phi (p, q, \hat{{\mathbf {p}}} \cdot \hat{{\mathbf {q}}}), \end{aligned}$$
(2.25)

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

$$\begin{aligned} \phi _{{\mathbf {p}}}({\mathbf {q}})= & {} \sum _l \frac{2l+1}{4\pi } \, F_l (p, q) P_l (\hat{{\mathbf {p}}} \cdot \hat{{\mathbf {q}}}) \nonumber \\= & {} \sum _{l, l_z} F_l (p, q) Y_{l}^{l_z \, *}(\theta _q, \phi _q) Y_l^{l_z}(\theta _p, \phi _p). \end{aligned}$$
(2.26)

As shown in Eq. (2.9), the perturbative amplitude can be expanded in spherical harmonics as well

$$\begin{aligned} {\mathcal {M}}^0_{\mathbf{m }_i \mathbf{m }_f}({\mathbf {q}})= & {} \sum _{l,l_z} \sum _{s, s_z} \langle s_1 m_1 s_2 m_2 | s s_z \rangle \nonumber \\&\times \,{\mathcal {M}}(q; l l_z; s s_z; \mathbf{m }_f) \, Y^{l_z}_l(\theta _q, \phi _q). \end{aligned}$$
(2.27)

Injecting Eqs. (2.26) and (2.27) into Eq. (2.23), the Sommerfeld-corrected matrix element can then be decomposed as

$$\begin{aligned} {\mathcal {M}}^{(S)}_{\mathbf{m }_i\mathbf{m }_f}({\mathbf {p}})= & {} \sum _{l, l_z} \sum _{s, s_z} \sum _{l', l_z'} \int \frac{{\mathrm {d}}^3 q}{(2\pi )^3} \, \langle s_1 m_1 s_2 m_2 | s s_z \rangle \nonumber \\&\times \,{\mathcal {M}}(q; l l_z; s s_z; \mathbf{m }_f) \, F_{l'} (p, q) \nonumber \\&\times \, Y^{l_z}_l(\theta _q, \phi _q) Y^{l_z'}_{l'}(\theta _q, \phi _q) Y^{l_z'}_{l'}(\theta _p, \phi _p) \nonumber \\= & {} \sum _{l, l_z} \sum _{s, s_z} \langle s_1 m_1 s_2 m_2 | s s_z \rangle \! \int \! \frac{q^2{\mathrm {d}} q}{2\pi ^2}\nonumber \\&\times \, {\mathcal {M}}(q; l l_z; s s_z; \mathbf{m }_f) \, F_l (p, q)\, Y^{l_z}_{l}(\theta _p, \phi _p).\nonumber \\ \end{aligned}$$
(2.28)

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:

$$\begin{aligned} {\mathcal {M}}^{(S)}_{ll_z;ss_z;\mathbf{m }_f}(p) = \int \frac{q^2{\mathrm {d}} q}{2\pi ^2}\, {\mathcal {M}}(q; l l_z; s s_z; \mathbf{m }_f) \, F_l (p, q).\nonumber \\ \end{aligned}$$
(2.29)

Using Eq. (2.12), we can re-express this amplitude as

$$\begin{aligned} {\mathcal {M}}^{(S)}_{ll_z;ss_z;\mathbf{m }_f}(p) = \sum _{n\ge 0} \alpha ^{\mathbf{m }_f}_{ll_z;ss_z;n} \int \frac{{\mathrm {d}} q}{2\pi ^2}\, q^{l + 2n + 2} \, F_l (p, q).\nonumber \\ \end{aligned}$$
(2.30)

As shown in [22], the integrals can be rewritten as functions of the derivatives of the radial components of the wave function \(R_{pl}(r)\),

$$\begin{aligned}&\int \frac{{\mathrm {d}} q}{2 \pi ^2}\, q^{l + 2n + 2} \, F_l (p, q)\nonumber \\&\quad =\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}. \end{aligned}$$
(2.31)

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

$$\begin{aligned} R_{pl}(z; x) = \mathrm{e}^\frac{\pi x}{4} \mathrm{e}^\frac{-iz}{2} z^l \sum _{j = 0}^\infty \frac{\Gamma \left( 1 + \frac{i x}{2} + l + j \right) }{(2l + 1 + j)!}\frac{(iz)^j}{j!},\nonumber \\ \end{aligned}$$
(2.32)

where \(z = 2 r p\), \(x = A m / 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 Eq. (2.32), we can then write

$$\begin{aligned} {\mathcal {M}}^{(S)}_{ll_z;ss_z;\mathbf{m }_f}(p) = \sum _{n\ge 0} \alpha ^{\mathbf{m }_f}_{ll_z;ss_z;n} p^{l + 2n} {\mathcal {C}}_{l} (x) {\mathcal {D}}_{ln} (x). \end{aligned}$$
(2.33)

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

$$\begin{aligned} {\mathcal {C}}_l (x)= & {} \frac{1}{(-i)^l} \, \mathrm{e}^\frac{\pi x}{4} \,\Gamma \left( 1 + \frac{i x}{2}\right) \prod _{b = 1}^{l} \left( 1 + \frac{i x}{2 b} \right) \nonumber \\ {\mathcal {D}}_{ln} (x)= & {} \frac{n! (2l+2n+1)!}{(l+n)!} \sum _{j = 0}^{2n} \frac{(-2)^j (l+j)!}{j! (2n-j)! (2l+j+1)!} \nonumber \\&\times \left[ \prod _{b = l+1}^{l+j}\left( 1 + \frac{i x}{2b}\right) \right] , \end{aligned}$$
(2.34)

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 Eq. (2.11) can then be written as

$$\begin{aligned} \left| {\mathcal {M}}^{(S)}_{ll_z;ss_z;\mathbf{m }_f}(p) \right| ^2= & {} S_l (x) \sum _{n,n'} \alpha ^{\mathbf{m }_f}_{ll_z;ss_z;n} \left( \alpha ^{\mathbf{m }_f}_{ll_z;ss_z;n'}\right) ^* \nonumber \\&\times \,{\mathcal {D}}_{ln} (x) {\mathcal {D}}_{ln'}^* (x) \, p^{2(l + n + n')}, \end{aligned}$$
(2.35)

where

$$\begin{aligned} S_{l} (x) = |{\mathcal {C}}_{l} (x)|^2 = \frac{\pi x}{1-\mathrm{e}^{-\pi x}} \prod _{b = 1}^{l} \left( 1 + \frac{x^2}{4 b^2} \right) \end{aligned}$$
(2.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)\) [21, 23]. Here, we used \(\left| \Gamma (1 + i b)\right| = \sqrt{\pi b \, \mathrm {csch}(\pi b)}\). Note that, in Eq. (2.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 Eq. (2.35) into Eq. (2.11).

2.3 Convergence and strategy

The Sommerfeld corrections as given in Eq. (2.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 Eq. (2.34), at low velocity, the Sommerfeld factor for a given (ln) is \({\mathcal {O}}(p^{-l - 2n - \frac{1}{2}})\). For a momentum expansion of the perturbative amplitude of the form

$$\begin{aligned} {\mathcal {M}}_{ll_z;ss_z;\mathbf{m }_f}^0(p) = \sum _{n\ge 0} \alpha _{ll_z;ss_z;n}^{\mathbf{m }_f} p^{l+2n}, \end{aligned}$$
(2.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 Eq. (2.33) becomes

$$\begin{aligned} {\mathcal {M}}_{ll_z;ss_z;\mathbf{m }_f}^{(S)}(p)= & {} \sqrt{\frac{\pi A m}{p}}\sum _{n\ge 0} (-1)^{l+n} \alpha _{ll_z;ss_z;n}^{\mathbf{m }_f} m^{l + 2n} \nonumber \\&\times \,\frac{A^{l+2n}}{2^l(l + n)!}\, \frac{n!}{(2n)!} + {\mathcal {O}}(p^\frac{1}{2}) \nonumber \\= & {} \sqrt{\frac{\pi Am}{p}}\sum _{n\ge 0} \widetilde{\alpha }_{ll_z;ss_z;n}^{\mathbf{m }_f}\nonumber \\&\times \,\frac{A^{l+2n}}{2^l(l + n)!}\, \frac{n!}{(2n)!} + {\mathcal {O}}(p^\frac{1}{2}), \end{aligned}$$
(2.38)

where \(\widetilde{\alpha }_{ll_z;ss_z;n}^{\mathbf{m }_f} \equiv (-1)^{l+n} \alpha _{ll_z;ss_z;n}^{\mathbf{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 (ln) 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 the Sommerfeld-corrected cross sections. The choice for \(l_{\mathrm {max}}\) determines the degree of precision for both expansions according to Eqs. (2.12) and (2.38).

  2. 2.

    For each value of l, include all expansion terms from Eq. (2.35) with \(n, n'\) satisfying \(n + n' + l \le 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 Eq. (2.35) into Eq. (2.11).

In this procedure the perturbative amplitude is fully expanded up to \(p^{l_{\mathrm {max}}}\) and the perturbative cross section up to \(p^{2 l_{\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 Sects. 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 Sect. 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 [16, 17] 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 [17, 34, 35]

$$\begin{aligned} V_{\mathrm {QCD}}= & {} C \frac{\alpha _s (\hat{\mu })}{r} \left[ 1 + \frac{\alpha _s (\hat{\mu }) }{4 \pi } \left( c_1 + 2 c_2 (\gamma _E + \log \hat{\mu } r) \right) \right] \nonumber \\\approx & {} C \frac{\alpha _s (\hat{\mu } \approx 1/r)}{r}, \end{aligned}$$
(3.1)

where C is proportional to the quadratic Casimir. For example for the quark–anti-quark 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 Eq. (3.1). Now, as shown in [16, 17], the QCD potential between two particles of SU(3) representations \(\mathbf{R }\) and \(\mathbf{R }'\) can be rewritten as a sum of Coulomb potentials of the form

$$\begin{aligned} V_{\mathbf{R } \otimes \mathbf{R }'}= & {} \frac{\alpha _s (\hat{\mu })}{r} \sum _a T^a_{\mathbf{R }} \otimes T^a_{\mathbf{R }'}\nonumber \\= & {} \frac{\alpha _s (\hat{\mu })}{2r} \sum _{\mathbf{Q }} \Big [ C_2(\mathbf{Q }) \mathbb {1}_{\mathbf{Q }} - C_2(\mathbf{R }) \mathbb {1} - C_2(\mathbf{R }') \mathbb {1} \Big ],\nonumber \\ \end{aligned}$$
(3.2)

where \(\mathbf{R } \otimes \mathbf{R }' = \bigoplus _{\mathbf{Q }} \mathbf{Q }\) and \(C_2(\mathbf{R })\), \(C_2(\mathbf{R }')\) are the quadratic Casimir indices for \(\mathbf{R }\) and \(\mathbf{R }'\), respectively. Each irreducible \(\mathbf{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 \(\mathbf{R } = \mathbf{3 }, \mathbf{6 }, \mathbf{8 }\) and \(\mathbf{R }'=\overline{\mathbf{R }}\). The corresponding color decompositions are

$$\begin{aligned} \begin{array}{ll} \mathbf{3 } \otimes \overline{\mathbf{3 }} &{} = \mathbf{1 } \oplus \mathbf{8 } \\ \mathbf{6 } \otimes \overline{\mathbf{6 }} &{} = \mathbf{1 } \oplus \mathbf{8 } \oplus \mathbf{27 } \\ \mathbf{8 } \otimes \mathbf{8 } &{} = \mathbf{1 }_\mathbf S \oplus \mathbf{8 }_\mathbf A \oplus \mathbf{8 }_\mathbf S \oplus \mathbf{10 }_\mathbf A \oplus \overline{\mathbf{10 }}_\mathbf A \oplus \mathbf{27 }_\mathbf S . \end{array} \end{aligned}$$
(3.3)

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

$$\begin{aligned} \begin{array}{c|c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \mathbf{R } &{} \mathbf{1 } &{} \mathbf{3 } &{} \mathbf{6 } &{} \mathbf{8 } &{} \mathbf{10 } &{} \mathbf{15 } &{} \mathbf{27 } &{} \mathbf{64 } \\ \hline C(\mathbf{R }) &{} 0 &{} \frac{1}{2} &{} \frac{5}{2} &{} 3 &{} \frac{15}{2} &{} 10 &{} 27 &{} 120 \\ C_2(\mathbf{R }) &{} 0 &{} \frac{4}{3} &{} \frac{10}{3} &{} 3 &{} 6 &{} \frac{16}{3} &{} 8 &{} 15 \\ \end{array} \end{aligned}$$
(3.4)

Injecting Eqs. (3.3) and (3.4) into Eq. (3.2), we find

Fig. 2
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

$$\begin{aligned} V_{\mathbf{3 } \otimes \overline{\mathbf{3 }}}= & {} \frac{\alpha _s}{r} \! \left\{ \!\! \begin{array}{ll} - \frac{4}{3} &{} \quad (\mathbf{1 }) \\ + \frac{1}{6} &{} \quad (\mathbf{8 }), \end{array} \right. \quad V_{\mathbf{6 } \otimes \overline{\mathbf{6 }}}=\frac{\alpha _s}{r} \! \left\{ \!\! \begin{array}{ll} - \frac{10}{3} &{} \quad (\mathbf{1 }) \\ - \frac{11}{6} &{} \quad (\mathbf{8 }) \\ + \frac{2}{3} &{} \quad (\mathbf{27 }), \end{array} \right. \nonumber \\ V_{\mathbf{8 } \otimes \mathbf{8 }}= & {} \frac{\alpha _s}{r} \! \left\{ \!\! \begin{array}{ll} - 3 &{} \quad (\mathbf{1 }_\mathbf S ) \\ - \frac{3}{2} &{} \quad (\mathbf{8 }_\mathbf A , \mathbf{8 }_\mathbf S ) \\ 0 &{} \quad (\mathbf{10 }_\mathbf A , \overline{\mathbf{10 }}_\mathbf A ) \\ + 1 &{} \quad (\mathbf{27 }_\mathbf S ). \end{array} \right. \end{aligned}$$
(3.5)

For a particle in a color representation \(\mathbf{R } = \mathbf{3 }, \mathbf{6 }, \mathbf{8 }\), the particle–antiparticle QCD potential at tree-level can be decomposed into Coulomb potentials with coupling strengths set by Eq. (3.5). The Coulomb interaction associated to a given irreducible representation \(\mathbf{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 (ls) expansion. To obtain the full Sommerfeld-corrected amplitude one has to find the irreducible representations \(\mathbf{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

$$\begin{aligned} \Phi \, \overline{\Phi } \rightarrow q_i \, \bar{q}_j \quad {\mathrm {and}} \quad \Phi \, \overline{\Phi } \rightarrow g^a \, g^b. \end{aligned}$$
(3.6)

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 Fig. 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 Appendices B.2 and B.3.

First we discuss the color structure of the amplitude for the annihilation into a quark–anti-quark pair. As shown in Fig. 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 \(\mathbf{R }\) of \(\Phi \):

$$\begin{aligned} {{\mathcal {A}}}^a \big |^i_j \propto (T_{\mathbf{R }}^a)^i_j, \end{aligned}$$
(3.7)

where a is the index of the s-channel gluon and the indices i and j run from 1 to the dimensionality of the \(\mathbf{R }\) representation, \(d_{\mathbf{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,

$$\begin{aligned} \sum _{{\mathrm {color}}} \left| A_{\mathbf{R } \otimes \overline{\mathbf{R }}} \right| ^2 = \sum _{{\mathrm {color}}} \big | [\mathbf{8 }] \big |^2, \end{aligned}$$
(3.8)

where \(\sum _{\mathrm {color}}\) runs over all the color indices of the external particles in the amplitude.Footnote 2 The octet representation \([\mathbf{8 }]\) is antisymmetric for the decomposition of self-conjugate representations, like the \(\mathbf{8 } \otimes \mathbf{8 }\), because of the CP nature of the exchanged gluon [36].

Annihilation processes into gluons have a more complex color structure. As can be seen in Fig. 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_{\mathbf{R }} \, T^b_{\mathbf{R }}\) and \(T^b_{\mathbf{R }} \, T^a_{\mathbf{R }}\) where ab are the color indices of the final-state gluons. In full generality, the amplitude can then be written as

$$\begin{aligned} {\mathcal {A}}^{ab} \big |^i_j = \alpha \left\{ T_{\mathbf{R }}^a, T_{\mathbf{R }}^b \right\} ^i_j + \beta \left[ T_{\mathbf{R }}^a, T_{\mathbf{R }}^b \right] ^i_j, \end{aligned}$$
(3.9)

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

$$\begin{aligned} {\mathcal {A}}^{ab} \big |^i_j = (-1)^{l+s} {\mathcal {A}}^{ba} \big |^i_j. \end{aligned}$$
(3.10)

The annihilation amplitude will therefore be proportional to the anticommutator of the \(T^a_{\mathbf{R }}\) for even \(l + s\) and to the commutator for odd \(l + s\).Footnote 3 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 Sect. 3.1, we consider particle–antiparticle annihilation with \(\mathbf{R } = \mathbf{3 }, \mathbf{6 }, \mathbf{8 }\). For amplitudes proportional to \([T^a_{\mathbf{R }}, T^b_{\mathbf{R }}]\), we can write

$$\begin{aligned} {\mathcal {A}}^{ab} \big |^i_j \propto [T^a_{\mathbf{R }}, T^b_{\mathbf{R }}]^i_j = i f^{abc} (T^c_{\mathbf{R }})^i_j. \end{aligned}$$
(3.11)

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

$$\begin{aligned} \sum _{\mathrm {color}} \left| A_{\mathbf{R } \otimes \overline{\mathbf{R }}} \right| ^2 = \sum _{\mathrm {color}} \big | [\mathbf{8 }] \big |^2. \end{aligned}$$
(3.12)

Here the \([\mathbf{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 Eq. (3.10).

For terms proportional to \(\left\{ T^a_{\mathbf{R }}, T^b_{\mathbf{R }} \right\} \), the amplitude decomposition depends on \(\mathbf{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 [38]. 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_{\mathbf{R }}, T^b_{\mathbf{R }} \right\} ^i_j\) tensor allows one to express the amplitude as

$$\begin{aligned} {\mathcal {A}}^{ab} \big |^i_j = \sum _{\mathbf{Q }} \, [\mathbf{Q }]^{ab} \big |^i_j \,, \end{aligned}$$
(3.13)

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

$$\begin{aligned} \sum _{\mathrm {color}} \big | A_{\mathbf{R } \otimes \overline{\mathbf{R }}} \big |^2 = \sum _{\mathbf{Q }} \left[ \sum _{\mathrm {color}} \big | [\mathbf{Q }] \big |^2 \right] . \end{aligned}$$
(3.14)

For the \(\Phi \, \overline{\Phi } \rightarrow g^a \, g^b\) process that we consider here, using Eqs. (B.1)–(B.3) for terms proportional to \(\left\{ T^a_{\mathbf{R }}, T^b_{\mathbf{R }} \right\} ^i_j\), we obtain the magnitude of the contributions from the different color states to the total amplitude. For \(\mathbf{R } = \mathbf{3 }, \mathbf{6 }, \mathbf{8 }\), we have

$$\begin{aligned} \sum _\mathrm{color} \big | A_{\mathbf{3 } \otimes {\bar{\mathbf{3}}}} \big |^2= & {} \frac{7}{2} \sum _\mathrm{color} \big | [\mathbf{1 }] \big |^2 = \frac{7}{5} \sum _\mathrm{color} \big | [\mathbf{8 }] \big |^2 \nonumber \\ \sum _\mathrm{color} \big | A_{\mathbf{6 } \otimes {\bar{\mathbf{6}}}} \big |^2= & {} \frac{31}{5} \sum _\mathrm{color} \big | [\mathbf{1 }] \big |^2 \nonumber \\= & {} \frac{155}{49} \sum _\mathrm{color} \big | [\mathbf{8 }] \big |^2 = \frac{155}{81} \sum _\mathrm{color} \big | [\mathbf{27 }] \big |^2 \nonumber \\ \sum _\mathrm{color} \big | A_{\mathbf{8 } \otimes \mathbf{8 }} \big |^2= & {} 6 \sum _\mathrm{color} \big | [\mathbf{1 }_{\mathbf{S}}] \big |^2 = 3 \sum _\mathrm{color} \big | [\mathbf{8 }_{\mathbf{S}}] \big |^2 \nonumber \\= & {} 2 \sum _\mathrm{color} \big | [\mathbf{27 }_{\mathbf{S}}] \big |^2. \end{aligned}$$
(3.15)

These results for the triplet and the octet agree with the ones obtained for the s-wave in [17, 39]. The results for the sextet and the more exotic decompositions discussed in Appendices B.2 and B.3 are novel and can also be used to extend the scope of the bound state calculations of [39] as described in [32].

3.3 Sommerfeld corrections

Combining the results from Sects. 3.1 and 3.2, the Sommerfeld-corrected cross sections for the annihilation of two colored states in the representations \(\mathbf{R }\) and \(\overline{\mathbf{R }}\) can be decomposed as

$$\begin{aligned} \sigma ^{(S)} = \sum _{\mathbf{Q }} \kappa _{\mathbf{Q }} \, \sigma ^{(S)}_{\mathrm {C}} \left[ \alpha _{\mathbf{Q }} \right] , \end{aligned}$$
(3.16)

where \(\mathbf{R } \otimes \overline{\mathbf{R }} = \bigoplus _{\mathbf{Q }} \mathbf{Q }\). \(\sigma ^{(S)}_C \left[ \alpha _{\mathbf{Q }} \right] \) is the Sommerfeld-corrected cross section for a Coulomb potential with coupling strength \(A = \alpha _{\mathbf{Q }}\), which can be computed by combining Eqs. (2.11) and (2.35). \(\kappa _{\mathbf{Q }}\) is the relative magnitude of the contribution of the \(\mathbf{Q }\) initial state to the annihilation amplitude, defined as

$$\begin{aligned} \sum _{\mathrm {color}} \big | [\mathbf{Q }] \big |^2 = \kappa _{\mathbf{Q }} \sum _{\mathrm {color}} \big | A_{\mathbf{R } \otimes \overline{\mathbf{R }}} \big |^2. \end{aligned}$$
(3.17)

As described in Sect. 3.2, the \(\kappa _{\mathbf{Q }}\) weights depend not only on the color representation of the initial state, but also on its (ls) 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 _{\mathbf{Q }}\) factors. In what follows, we will therefore consider cross sections associated to an individual (ls) particle–antiparticle initial state in the \(\mathbf{R } \otimes \overline{\mathbf{R }}\) representation.

Reading off \(\alpha _{\mathbf{Q }}\) from Eq. (3.5) and \(\kappa _{\mathbf{Q }}\) from Eqs. (3.8), (3.12) and (3.15), for \(\mathbf{R } = \mathbf{3 }, \mathbf{6 }, \mathbf{8 }\), the Sommerfeld-corrected cross sections are

$$\begin{aligned} \sigma ^{(S)}_{\mathbf{3 } \otimes \overline{\mathbf{3 }} \rightarrow q \, \overline{q}}= & {} \sigma ^{(S)}_C \left[ - \frac{\alpha _s}{6} \right] , \nonumber \\ \sigma ^{(S)}_{\mathbf{3 } \otimes \overline{\mathbf{3 }} \rightarrow g \, g}= & {} {\left\{ \begin{array}{ll} \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 ,\nonumber \\ \sigma ^{(S)}_C \left[ - \frac{\alpha _s}{6} \right] &{} {\mathrm {odd}} \,\, l + s ,\end{array}\right. } \nonumber \\ \sigma ^{(S)}_{\mathbf{6 } \otimes \overline{\mathbf{6 }} \rightarrow q \, \overline{q}}= & {} \sigma ^{(S)}_C \left[ \frac{11 \alpha _s}{6} \right] ,\nonumber \\ \sigma ^{(S)}_{\mathbf{6 } \otimes \overline{\mathbf{6 }} \rightarrow g \, g}= & {} {\left\{ \begin{array}{ll} \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] \nonumber \\ \quad + \frac{81}{155} \sigma ^{(S)}_C \left[ - \frac{2 \alpha _s}{3} \right] &{} {\mathrm {even}} \,\, l + s ,\nonumber \\ \sigma ^{(S)}_C \left[ \frac{11 \alpha _s}{6} \right] &{} {\mathrm {odd}} \,\, l + s, \end{array}\right. } \nonumber \\ \sigma ^{(S)}_{\mathbf{8 } \otimes \mathbf{8 } \rightarrow q \, \overline{q}}= & {} \sigma ^{(S)}_C \left[ \frac{3 \alpha _s}{2} \right] ,\nonumber \\ \sigma ^{(S)}_{\mathbf{8 } \otimes \mathbf{8 } \rightarrow g \, g}= & {} {\left\{ \begin{array}{ll} \frac{1}{6} \sigma ^{(S)}_C \left[ 3 \alpha _s \right] {+} \frac{1}{3} \sigma ^{(S)}_C \left[ \frac{3 \alpha _s}{2} \right] \\ \quad + \frac{1}{2} \sigma ^{(S)}_C\left[ - \alpha _s \right] &{} \mathrm {even} \,\, l + s \\ \sigma ^{(S)}_C \left[ \frac{3 \alpha _s}{2} \right] &{} {\mathrm {odd}} \,\, l + s .\end{array}\right. }.\nonumber \\ \end{aligned}$$
(3.18)

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 Sect. 2. The final analytic expressions for the Sommerfeld-corrected cross sections can be found by combining Eqs. (2.11) and (2.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 [40, 41], 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 Sect. 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 Eq. (3.5) for the QCD potential. This new potential can be straightforwardly obtained from Eq. (3.18) by averaging over the different channels. We then obtain

$$\begin{aligned}&V_{\mathbf{3 } \otimes \overline{\mathbf{3 }} \rightarrow g \, g}^{{\mathrm {avg}}} = - \frac{11}{42} \frac{\alpha _s}{r}, \quad V_{\mathbf{6 } \otimes \overline{\mathbf{6 }} \rightarrow g \, g}^{{\mathrm {avg}}} = - \frac{143}{186} \frac{\alpha _s}{r}, \nonumber \\&\quad V_{\mathbf{8 } \otimes \mathbf{8 } \rightarrow g \, g}^{{\mathrm {avg}}} = - \frac{1}{2} \frac{\alpha _s}{r}. \end{aligned}$$
(3.19)

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

$$\begin{aligned} \sigma ^{(S), \, {\mathrm {avg}}}_{\mathbf{3 } \otimes \overline{\mathbf{3 }} \rightarrow g \, g}= & {} \sigma ^{(S)}_C \left[ \frac{11 \alpha _s}{42} \right] \nonumber \\ \sigma ^{(S), \, {\mathrm {avg}}}_{\mathbf{6 } \otimes \overline{\mathbf{6 }} \rightarrow g \, g}= & {} \sigma ^{(S)}_C \left[ \frac{143 \alpha _s}{186} \right] \\ \sigma ^{(S), \, {\mathrm {avg}}}_{\mathbf{8 } \otimes \mathbf{8 } \rightarrow g \, g}= & {} \sigma ^{(S)}_C \left[ \frac{\alpha _s}{2} \right] .\nonumber \end{aligned}$$
(3.20)

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 [32], 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

(4.1)

where ij are color indices and the \(T^a_{\mathbf{R }}\) matrices are the generators for the color representation \(\mathbf{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

$$\begin{aligned} \begin{aligned} V^{\mu \nu }_i&= D^\mu _{ij} V^\nu _j - D^\nu _{ij} V^\mu _j , \\ D_{\mu , ij}&= \partial _\mu \delta _{ij} - i g_s G_\mu ^a (T^a_{\mathbf{R }})_{ij}. \end{aligned} \end{aligned}$$
(4.2)

Note that the Lagrangian for vectors can also include anomalous terms [42, 43] 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 [32].

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

$$\begin{aligned}&\sigma \! (S \, S \! \rightarrow \! q \, \bar{q})=\frac{2 \pi \alpha _s^2}{3 s} \frac{C_2(\mathbf{R })}{d_{\mathbf{R }}} \beta _S ,\nonumber \\&\sigma \! (S \, S \! \rightarrow \! g \, g)=\frac{2 \pi \alpha _s^2}{3 s^3} \frac{C_2(\mathbf{R })}{d_{\mathbf{R }} \beta _S^2} \nonumber \\&\quad \times \left[ C_2(\mathbf{G })\left( s \beta _S (10 m_S^2 - s) - 24 m_S^4 \log \frac{1+\beta _S}{1-\beta _S}\right) \right. \nonumber \\&\quad \left. +\, 6 C_2(\mathbf{R })\left( s \beta _S (s {+} 4 m_S^2) {+} m_S^2 (8 m_S^2 - 4 s) \log \frac{1 \! + \! \beta _S}{1 \! - \! \beta _S} \right) \! \right] ,\nonumber \\&\sigma (\psi \, \bar{\psi } \! \rightarrow \! q \, \bar{q})= \frac{2 \pi \alpha _s^2}{3 s} \frac{C_2(\mathbf{R })}{d_{\mathbf{R }}} \frac{1}{\beta _\psi } \left( 1 + \frac{2 m_\psi ^2}{s}\right) , \nonumber \\&\sigma (\psi \, \bar{\psi } \! \rightarrow \! g \, g)= - \frac{2 \pi \alpha _s^2}{3 s^3} \frac{C_2(\mathbf{R })}{d_{\mathbf{R }} \beta _\psi ^2} \nonumber \\&\quad \times \left[ C_2(\mathbf{G }) \left( s \beta _\psi (s + 5 m_\psi ^2) - 12 m_\psi ^4 \log \frac{1 \! + \! \beta _\psi }{1 \! - \! \beta _\psi } \right) \right. \nonumber \\&\quad +\, 3 C_2(\mathbf{R }) \left( s \beta _\psi (s + 4 m_\psi ^2) \phantom {\frac{1}{2}}\right. \nonumber \\&\quad \left. \left. +\,(8 m_\psi ^4 - 4 m_\psi ^2 s - s^2) \log \frac{1 \! + \! \beta _\psi }{1 \! - \! \beta _\psi } \right) \right] ,\nonumber \\&\sigma ( V \, V \! \rightarrow \! q \, \bar{q})=\frac{\pi \alpha _s^2}{54 s} \frac{C_2(\mathbf{R })}{d_{\mathbf{R }}} \beta _V \frac{12 m_V^4 + 20 m_V^2 s + s^2}{m_V^4} \nonumber \\&\sigma (V \, V \! \rightarrow \! g \, g)=\frac{2 \pi \alpha _s^2}{9 m_V^2 s^3} \frac{C_2(\mathbf{R })}{d_{\mathbf{R }}\beta _V^2} \nonumber \\&\quad \times \left[ C_2(\mathbf{G }) m_V^2 \!\! \left( \! s \beta _V \! (10 m_V^2 \! + \! 7 s ) \! - \! 8 \! ( 3 m_V^4 \! + \! s^2) \log \frac{1 \! + \! \beta _V}{1 \! - \! \beta _V} \! \right) \! \right. \nonumber \\&\quad +\,2 C_2(\mathbf{R }) \left( s \beta _V \! ( 12 m_V^4 \! + \! 3 m_V^2 s \! + \! 4 s^2 )\right. \nonumber \\&\quad \left. \left. +\, 12 ( 2 m_V^6 \! - \! m_V^4 s ) \log \frac{1 \! + \! \beta _V}{1 \! - \! \beta _V} \right) \right] . \end{aligned}$$
(4.3)

In these expressions, the phase space factor is defined by \(\beta _\Phi = \sqrt{1 - \frac{4 m_\Phi ^2}{s}}\) and \(C_2(\mathbf{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 [32].

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 Sects. 2 and 3 as well as the perturbative amplitudes for the models described in Sect. 4 in a Mathematica notebook that is attached to this paper [44]. This notebook also provides an interface to micrOMEGAs [45, 46] 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.

Fig. 3
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}} \le m_\Phi \le 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 [47, 48].

The results for different annihilation processes are shown in Fig. 3. To outline the mass dependence of the Sommerfeld-corrected ratios discussed before, we plot these ratios as a band for \(500~\mathrm {GeV} \le m_\Phi \le 2500~\mathrm {GeV}\). We first notice that, as mentioned in Sect. 2, in spite of the \({\mathcal {O}}(v^{-1})\) terms present at large l due to Sommerfeld corrections both the perturbative and the 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 Fig. 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 [32], 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 Sect. 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 [24], hidden valleys [49, 50], dark radiation [25, 26] and glueball dark matter [27,28,29,30,31] 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 Sects. 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 [27,28,29,30,31]. 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 Sects. 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 [31] 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 Sect. 3 to particles charged under a new dark gauge group SU(N), either in the fundamental \(\mathbf{F }\) or in the adjoint \(\mathbf{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)

$$\begin{aligned} \Phi \, \overline{\Phi } \rightarrow Q_i \, \bar{Q}_j \quad {\mathrm {and}} \quad \Phi \, \overline{\Phi } \rightarrow G^a \, G^b. \end{aligned}$$
(5.1)

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 Sect. 3. In particular, the leading-order term of the SU(N) potential can be described by a Coulomb potential obeying Eq. (3.2) and the symmetry constraints on the different partial waves arising from CP conservation are independent on N.

A generalized version of Eq. (3.3) can be derived by decomposing the \(\mathbf{F } \otimes \overline{\mathbf{F }} \) and the \(\mathbf{A } \otimes \mathbf{A }\) products to obtain the following possible representations for the \(\Phi \, \overline{\Phi }\) initial state:

$$\begin{aligned} \begin{array}{ll} \mathbf{F } \otimes \overline{\mathbf{F }} &{} = \mathbf{1 } \oplus \mathbf{A } \\ \mathbf{A } \otimes \mathbf{A } &{} = \mathbf{1 }_\mathbf S \oplus \mathbf{A }_\mathbf A \oplus \mathbf{A }_\mathbf S \oplus \mathbf{B }_\mathbf S \oplus \mathbf{C }_\mathbf A \oplus \overline{\mathbf{C }}_\mathbf A \oplus \mathbf{D }_\mathbf S \,. \end{array} \end{aligned}$$
(5.2)

A notable difference from the SU(3) case here is the appearance of the \(\mathbf{B }_\mathbf S \) representation for \(N\ge 4\). The representations in Eq. (5.2) are associated with the following Young tableaux:

(5.3)

where \(\mathbf{A }\) and \(\mathbf{D }\) have \(N-1\) vertical boxes and \(\mathbf{B }\) and \(\mathbf{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 [51, 52] are summarized:

$$\begin{aligned} \begin{array}{c | c c c c c c} \mathbf{R } &{} \mathbf{1 } &{} \mathbf{F } &{} \mathbf{A } &{} \mathbf{B } &{} \mathbf{C } &{} \mathbf{D } \\ \hline \mathrm {dim}(\mathbf{R }) &{} 1 &{} N &{} N^2 \! - \! 1 &{} \frac{1}{4} (N^4 \! - \! 2 N^3 \! - \! 3 N^2) &{} \frac{1}{4} (N^4 \! - \! 5 N^2 \! + \! 4) &{} \frac{1}{4} (N^4 \! + \! 2 N^3 \! - \! 3 N^2) \\ C(\mathbf{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(\mathbf{R }) &{} 0 &{} \frac{N^2 - 1}{2N} &{} N &{} 2 (N \! - \! 1) &{} 2 N &{} 2 (N \! + \! 1) \\ \end{array} \end{aligned}$$
(5.4)

Calling the new gauge coupling \(\alpha _N\), we can now use Eqs. (3.2) and (5.2) as well as the table in Eq. (5.4) to derive the SU(N) Coulomb potential associated with the different \(\Phi \, \overline{\Phi }\) representations:

$$\begin{aligned}&V_{\mathbf{F } \otimes \overline{\mathbf{F }}} = \! \frac{\alpha _N}{r} \! \left\{ \!\! \begin{array}{ll} - \frac{N^2 - 1}{2N} &{} \quad (\mathbf{1 }) \\ \frac{1}{2N} &{} \quad (\mathbf{A }) \end{array} \right. ,\nonumber \\&V_{\mathbf{A } \otimes \mathbf{A }} = \! \frac{\alpha _N}{r} \! \left\{ \!\! \begin{array}{ll} - N &{} \quad (\mathbf{1 }_\mathbf S ) \\ - \frac{N}{2} &{} \quad (\mathbf{A }_\mathbf A , \mathbf{A }_\mathbf S ) \\ - 1 &{} \quad (\mathbf{B }_\mathbf S ) \\ 0 &{} \quad (\mathbf{C }_\mathbf A , \overline{\mathbf{C }}_\mathbf A ) \\ 1 &{} \quad (\mathbf{D }_\mathbf S ). \end{array} \right. \end{aligned}$$
(5.5)

For the case of \(N = 3\), this potential reduces to Eq. (3.5). 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 Sect. 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 } \rightarrow 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 } \rightarrow G^a \, G^b\) annihilation process, the CP conservation arguments described in Sect. 3.2 still apply and, as in Eq. (3.12), the squared amplitude for odd \(l + s\) can be written as

$$\begin{aligned} \sum _{\mathrm {color}} \left| A_{\mathbf{R } \otimes \overline{\mathbf{R }}} \right| ^2 = \sum _{\mathrm {color}} \big | [\mathbf{A }] \big |^2, \end{aligned}$$
(5.6)

for all \(\mathbf{R }\). Similarly, for even \(l + s\), the decompositions given in Eq. (3.15) for the products of two fundamentals and two adjoints can be generalized to

$$\begin{aligned} \sum _{\mathrm {color}} \big | A_{\mathbf{F } \otimes \overline{\mathbf{F }}} \big |^2 = \frac{N^2 - 2}{2} \sum _{\mathrm {color}} \big | [\mathbf{1 }] \big |^2 = \frac{N^2 - 2}{N^2 - 4} \sum _{\mathrm {color}} \big | [\mathbf{A }] \big |^2 \nonumber \\ \end{aligned}$$
(5.7)

and

$$\begin{aligned} \begin{aligned} \sum _{\mathrm {color}} \big | A_{\mathbf{A } \otimes \mathbf{A }} \big |^2&= \frac{3}{4} (N^2 - 1) \sum _{\mathrm {color}} \big | [\mathbf{1 }_\mathbf S ] \big |^2\\ \sum _{\mathrm {color}} \big | A_{\mathbf{A } \otimes \mathbf{A }} \big |^2&= 3 \sum _{\mathrm {color}} \big | [\mathbf{A }_\mathbf S ] \big |^2 \\ \sum _{\mathrm {color}} \big | A_{\mathbf{A } \otimes \mathbf{A }} \big |^2&= \frac{3(N - 1)}{N-3} \sum _{\mathrm {color}} \big | [\mathbf{B }_\mathbf S ] \big |^2\\ \sum _{\mathrm {color}} \big | A_{\mathbf{A } \otimes \mathbf{A }} \big |^2&= \frac{3(N + 1)}{N+3} \sum _{\mathrm {color}} \big | [\mathbf{D }_\mathbf S ] \big |^2. \end{aligned} \end{aligned}$$
(5.8)

Note that these results only apply for \(N \ge 4\). For \(N = 3\) the contribution from the \(\mathbf{B }_\mathbf S \) representation goes to zero. In the case of \(N = 2\), if \(\Phi \) is in the fundamental representation the \(\Phi \, \overline{\Phi } \rightarrow 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 \(\mathbf{1 }_\mathbf S \) and the \(\mathbf{D }_\mathbf S = \mathbf{5 }_\mathbf S \) representation will contribute to the \(\Phi \, \overline{\Phi } \rightarrow 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 \(\mathbf{A }_\mathbf S \), \(\mathbf{B }_\mathbf S \) and \(\mathbf{D }_\mathbf S \) channels for initial-state particles in the adjoint representation.

5.2 Sommerfeld corrections

We now use the results from Sect. 5.1 as well as the methodology described in Sect. 3.3 to derive the Sommerfeld correction factors for the \(\Phi \, \overline{\Phi } \rightarrow Q_i \, \bar{Q}_j\) and the \(\Phi \, \overline{\Phi } \rightarrow G^a \, G^b\) annihilation processes. For general N, these factors now read

$$\begin{aligned} \sigma ^{(S)}_{\mathbf{F } \otimes \overline{\mathbf{F }} \rightarrow Q\, \overline{Q}}= & {} \sigma ^{(S)}_C \left[ - \frac{\alpha _N}{2N} \right] \nonumber \\ \sigma ^{(S)}_{\mathbf{F } \otimes \overline{\mathbf{F }} \rightarrow G \, G}= & {} {\left\{ \begin{array}{ll} \frac{2}{N^2 - 2} \sigma ^{(S)}_C \left[ \frac{(N^2 - 1) \alpha _N}{2N} \right] \\ \qquad + \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{array}\right. }\nonumber \\ \sigma ^{(S)}_{\mathbf{A } \otimes \mathbf{A } \rightarrow Q \, \overline{Q}}= & {} \sigma ^{(S)}_C \left[ \frac{N \alpha _N}{2} \right] ,\nonumber \\ \sigma ^{(S)}_{\mathbf{A } \otimes \mathbf{A } \rightarrow G \, G}= & {} {\left\{ \begin{array}{ll} \frac{4}{3 (N^2 - 1)} \sigma ^{(S)}_C \left[ N \alpha _N \right] \\ \qquad + \frac{1}{3} \sigma ^{(S)}_C \left[ \frac{N \alpha _N}{2} \right] &{} \\ \qquad + \frac{N-3}{3(N - 1)} \sigma ^{(S)}_C\left[ \alpha _N \right] \\ \qquad + \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{array}\right. } \end{aligned}$$
(5.9)

The ratios of the s, p and d-wave annihilation cross sections are shown in Fig. 4 for the \(\Phi \, \overline{\Phi } \rightarrow Q_i \, \bar{Q}_j\) and \(\Phi \, \overline{\Phi } \rightarrow 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 Fig. 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 \le N \le 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.

Fig. 4
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 \le N \le 10\) for a specific process and representation (either \(\mathbf{F }\) or \(\mathbf{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 Eq. (5.9), 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 Fig. 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 [31, 49, 50]. 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 [53]. The total potential is given as

$$\begin{aligned} V = \frac{\alpha _s (\hat{\mu })}{r} \sum _a T^a_{\mathbf{R }} \otimes T^a_{\mathbf{R }'} + \frac{\alpha _N (\hat{\mu })}{r} \sum _a T^a_{\mathbf{P }} \otimes T^a_{\mathbf{P }'}, \end{aligned}$$
(5.10)

where \(\mathbf{R }\) is a representation of SU(3) and \(\mathbf{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 Eq. (3.2) to each of the terms on the right-hand side of Eq. (5.10) 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 Eqs. (5.6)–(5.8). 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.

Fig. 5
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 Sect. 3.3. As outlined in Sect. 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

$$\begin{aligned} V = V_{SU(3)} - \frac{N^2 - 1}{2 N} \frac{\alpha _N}{r}, \end{aligned}$$
(5.11)

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

$$\begin{aligned} V = - \frac{4}{3} \frac{\alpha _s}{r} + V_{SU(N)}, \end{aligned}$$
(5.12)

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

Finally, the \(\psi \, \bar{\psi } \rightarrow g \, G\) process has not been studied before and has been not been taken into account in [31]. 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

$$\begin{aligned} V = \frac{1}{6} \frac{\alpha _s}{r} + \frac{1}{2 N} \frac{\alpha _N}{r}. \end{aligned}$$
(5.13)

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

$$\begin{aligned} (\sigma v)_{\psi \, \bar{\psi } \rightarrow q \, \bar{q}}= & {} 6 \times \frac{\pi \alpha _s^2}{9 N m_Q^2} \times S \left( - \frac{\alpha _s}{6 \beta } + \frac{N^2 \! - \! 1}{2 N} \frac{\alpha _N}{\beta } \right) , \nonumber \\ (\sigma v)_{\psi \, \bar{\psi } \rightarrow g \, g}= & {} \frac{7 \pi \alpha _s^2}{54 N m_Q^2} \times \left[ \frac{2}{7} S \left( \frac{4 \alpha _s}{3 \beta } + \frac{N^2 \! - \! 1}{2 N} \frac{\alpha _N}{\beta }\right) \right. \nonumber \\&\left. +\,\frac{5}{7} S \left( -\frac{\alpha _s}{6 \beta } + \frac{N^2 \! - \! 1}{2 N}\frac{\alpha _N}{\beta }\right) \right] ,\nonumber \\ (\sigma v)_{\psi \, \bar{\psi } \rightarrow G \, G}= & {} \frac{(N^2 \! - \! 1)(N^2 \! - \! 2) \pi \alpha _N^2}{48 N^3 m_Q^2} \nonumber \\&\times \left[ \frac{2}{N^2 \! - \! 2} S \left( \frac{4 \alpha _s}{3 \beta } + \frac{N^2 \! - \! 1}{2 N} \frac{\alpha _N}{\beta } \right) \right. \nonumber \\&\left. +\, \frac{N^2 \! - \! 4}{N^2 \! - \! 2} S \left( \frac{4 \alpha _s}{3 \beta } - \frac{1}{2N} \frac{\alpha _N}{\beta } \right) \right] ,\nonumber \\ (\sigma v)_{\psi \, \bar{\psi } \rightarrow g \, G}= & {} \frac{2 (N^2 \! - \! 1) \pi \alpha _s \alpha _N}{9 N^2 m_Q^2} \times S \left( {-} \frac{\alpha _s}{6 \beta } {-} \frac{1}{2N} \frac{\alpha _N}{\beta }\right) ,\nonumber \\ \end{aligned}$$
(5.14)

where the Sommerfeld factor S(x) is given in Eq. (2.4). The ratios of these cross sections over the s-wave perturbative cross sections for each process are shown in Fig. 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 2 to 8 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 } \rightarrow 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 [27,28,29,30,31].

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 (ls) 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 (ls) 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 Sommerfeld-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 (ls) 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 easily be included using the presented formalism.