A Sommerfeld Toolbox for Colored Dark Sectors

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.


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 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 swave 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 nonperturbative 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 , which can be either a scalar, a fermion or a vector. We first consider scenarios where 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 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 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.

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 → 2 process with two scalar fields in the initial state, this expansion would be of the form where p is the magnitude of the initial-state momentum in the center-of-mass frame, (θ, φ) are the scattering angles and Y lm (θ, φ) 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 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, 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] S where β is the velocity of the incoming particles in the centerof-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 ≥ 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 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.

Partial wave expansion
For a given field , the → SM 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 Sect. 2.2, obtaining Sommerfeld corrections involves computing the non-relativistic wave function for the two 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 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 (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 → g a g b in particular, since the gluons are identical particles, the contributions from states with even l + s will be proportional to T a R , T b R while the ones for states with odd l +s will be proportional to T a where T a R is the generator for the representation R of . 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 → 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ẑ-axis and denote the scattering angles in the initial state by (θ, φ). With {m 1 , m 2 } and {m 3 , m 4 } being the individual spin projections on theẑ-axis in the initial and final states, respectively, the total annihilation amplitude is defined as an element of the transition matrix T : (2.6) A given | p i ; ll z ; m 1 m 2 state can be decomposed into | p; ll z ; ss z states using Clebsch-Gordan coefficients, where s 1 , s 2 are the total spins of the incoming particles. The total amplitude can then be written as The matrix element 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 to the total amplitude. Knowing T f i , this matrix element can be computed using Since the (l, s) components of the amplitude are orthogonal the total cross section is of the form where d is the number of degrees of freedom of the field and d R is the dimensionality of the color representation of . Another factor 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 = s 4 − m 2 , with the lowest-order contribution for a given l being O( p l ). We can therefore write ll z ss z ,n p l+2n . (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).

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 → n process with the Sommerfeld effect occurring in the initial state, the amplitude then verifies the following recursion relation [33]: where α and β are the initial and final states respectively and the integral over γ represents the sum over all possible intermediate states. M 0 βα is the perturbative scattering amplitude corresponding to the exchange of one particle and V γ α is the non-relativistic interaction potential distorting the initial state α.
The interaction potential V γ α can be rewritten as where {m 1 , m 2 } and {m a , m b } are the z-components of the spins of the α and γ states respectively and p, k = |p|, |k| are the magnitudes of the momenta p and k in these states.
In the rest of this work, we will focus on a spin-independent spherically symmetric potential V (|p−k|). We can therefore factor out the spin states, which gives For initial-and final-state spins m i = {m 1 , m 2 } and m f = {m 3 , m 4 }, the Sommerfeld-corrected amplitude can be expressed as In the non-relativistic limit, the denominator can be rewritten as where E is the total energy of the system and μ is its reduced mass. For two particles of identical mass m, μ = m 2 . If the initial states are off-shell, that is, p 2 2μ = E, we can define [23] which verifies In position space (we go from p to r) this becomes where The final amplitude can now be computed by putting the initial states back on-shell, which leads to [21,23] where φ p (q) obeys the traditional Schrödinger equation in position space, For a potential of the form V (|p − q|), the wave function can be rewritten as φ p (q) = φ( p, q,p ·q), (2.25) and can therefore be expanded in Legendre polynomials and in spherical harmonics: As shown in Eq. (2.9), the perturbative amplitude can be expanded in spherical harmonics as well 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: 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 ), dq 2π 2 q l+2n+2 F l ( p, q) For a Coulomb potential V = −A/r , the radial components of the wave function can be computed analytically and are equal to where z = 2r p, 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 α s or the SU (N ) coupling α N . Using the expression given in Eq. (2.32), we can then write The Sommerfeld factors C l (x) and D ln (x) are given by where C l (x) is the correction to the amplitude for a perturbative matrix element of the form p l Y l z l (θ p , φ p ). Note here that 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 is the Sommerfeld correction for a Coulomb potential and for a perturbative amplitude of the form p l Y l z l (θ p , φ p ) [21,23].
. 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).

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 O(v l ) or equivalently 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 (l, where α m f ll z ;ss z ;n ≡ (−1) l+n α m f ll z ;ss z ;n 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 − 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: In this procedure the perturbative amplitude is fully expanded up to p l max and the perturbative cross section up to p 2l max −1 . Applying Sommerfeld corrections to this expansion gives an angular momentum expansion of the final cross section up to l 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.

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.

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] where C is proportional to the quadratic Casimir. For example for the quark-anti-quark potential C = 4 3 and the one-loop coefficients are defined by c 1 = 31 3 − 10 9 n f and c 2 = 11 − 2 3 n f , where n f is the number of active quark flavors at the scaleμ. it shows that the QCD potential at higher orders can be approximated as a simple Coulomblike 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 R and R can be rewritten as a sum of Coulomb potentials of the form where R ⊗ R = Q Q and C 2 (R), C 2 (R ) are the quadratic Casimir indices for R and R , 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 α s (μ) 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 α s in the rest of this section.
In what follows, we will consider particle-antiparticle annihilation with R = 3, 6, 8 and R = R. The corresponding color decompositions are The subscripts S and A indicate whether the representation is symmetric or antisymmetric, respectively, under the interchange of the two equal representations R and R . The quadratic Casimir indices (C 2 ) of these representations are given in Eq. (1) , For a particle in a color representation R = 3, 6, 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 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 Fig. 2 Feynman diagrams for the annihilation of into either a quark-anti-quark pair or a pair of gluons. The annihilating field may be scalar, fermion or vector, however, in the case of the fermion the four-point interaction is absent 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.

Decomposing perturbative cross sections
In this section, we consider tree-level annihilation of a particle into quarks and gluons through the strong interaction → q iq j and → g a g b . (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 . Here, we take 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 R of : 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 R . Only the color-octet configuration of the initial state, matching the representation of the exchanged gluon, will therefore contribute to the → q iq j cross section, where color runs over all the color indices of the external particles in the amplitude. 2 The octet representation [8] is antisymmetric for the decomposition of self-conjugate representations, like the 8 ⊗ 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 combi- are the color indices of the final-state gluons. In full generality, the amplitude can then be written as where α, β are factors that contain the kinematic dependence. As underlined in Sect. 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 The annihilation amplitude will therefore be proportional to the anticommutator of the T a R for even l+s and to the commutator for odd l + s. 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 → 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 R = 3, 6, 8. For amplitudes proportional to As for annihilation into qq, the amplitudes here are proportional to linear combinations of the generators of the R representation and therefore receive contributions from color-octet configurations only Here the [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 T a R , T b R , 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 [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 T a R , T b R i j tensor allows one to express the amplitude as for R⊗R = Q Q and where [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 (3.14) For the → g a g b process that we consider here, using 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].

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 R and R can be decomposed as α Q is the Sommerfeldcorrected cross section for a Coulomb potential with coupling strength A = α Q , which can be computed by combining Eqs. (2.11) and (2.35). κ Q is the relative magnitude of the contribution of the Q initial state to the annihilation amplitude, defined as color (3.17) As described in Sect. 3.2, the κ 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 → 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 κ Q factors.
In what follows, we will therefore consider cross sections associated to an individual (l, s) particle-antiparticle initial state in the R ⊗ R representation.
Reading off α Q from Eq.  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. 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.

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 with arbitrary spin-scalar, fermion or vector-and with an arbitrary representation under SU (3). We then introduce a small coupling between DM and 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 . 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 and compute the associated Sommerfeld corrections.

Simplified models
We consider scenarios where is either a real or complex scalar, Dirac or Majorana fermion or a real or complex vector boson. The kinetic and mass terms for = {S, ψ, V } in the complex scalar, Dirac fermion and complex vector models are then where i, j are color indices and the T a R matrices are the generators for the color representation R of . 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 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 , especially on perturbative unitarity, are discussed in our companion paper [32].
We list here the analytic cross sections for the pair annihilation to qq and g g. The total annihilation cross sections for = S, ψ, V are In these expressions, the phase space factor is defined by s and C 2 (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 in the Lagrangian, the V V → qq cross section grows as 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].

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 freezeout, 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.
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 α s (μ) 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 α Sommerfeld s (μ) 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 α 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 Sommerfeldcorrected ratios discussed before, we plot these ratios as a band for 500 GeV ≤ m ≤ 2500 GeV. We first notice that, as mentioned in Sect. 2, in spite of the 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 dwave perturbative cross section is indistinguishable from the exact value up to v ∼ 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 O(10%) of the Sommerfeldcorrected cross section. Although in several models these effects can be mitigated by a cancellation between the qq 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 sce- the Sommerfeld-corrected ratios are shown as a band corresponding to 500 GeV ≤ m ≤ 2500 GeV. For each of the processes we show the results for a specific color representation, which is denoted by the subscript on the fields narios where dark matter coannihilates with a colored dark partner.

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 α 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 ).

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 F or in the adjoint A representation. As before, we consider the self-annihilation of a particle into two fermions in the fundamental representation of SU (N ) or into two dark gauge bosons in the adjoint representation of SU (N ) 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 F ⊗ F and the A ⊗ A products to obtain the following possible representations for the initial state: A notable difference from the SU (3) case here is the appearance of the B S representation for N ≥ 4. The representations in Eq. (5.2) are associated with the following Young tableaux: . . . . . , (5.3) 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 [51,52] are summarized: 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 → Q iQ 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 → 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 for all 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 (5.8) Note that these results only apply for N ≥ 4. For N = 3 the contribution from the B S representation goes to zero. In the case of N = 2, if is in the fundamental representation the → G a G b process occurs only when is an SU (2) singlet. When is in the adjoint representation, only the states in the 1 S and the D S = 5 S representation will contribute to the → 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 A S , B S and D S channels for initial-state particles in the adjoint representation.

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 → Q iQ j and the → G a G b annihilation processes. For general N , these factors now read The ratios of the s, p and d-wave annihilation cross sections are shown in Fig. 4 for the → Q iQ j and → G a G b processes with 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 dwave but this time, we show the values of these different cross sections for 4 ≤ N ≤ 10. Contrary to the SU (3) case, we do not evaluate α N at the scale of the momenta of the incoming particles and instead set the coupling entering into the Sommerfeld corrections to be α Sommerfeld N (μ) = 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.
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 σ 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 σ C N α N 2 . 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.

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 Stan- Fig. 4 Ratios of the perturbative (solid lines) and Sommerfeldcorrected 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 ≤ N ≤ 10 for a specific process and representation (either F or A), which is denoted by the subscript on the fields dard 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 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 ( 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 ).

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.
In what follows, we consider a fermionic messenger particle ψ charged as a triplet under QCD and as a fundamental under SU (N ). In this scenario, ψ can annihilate either to g g, qq, 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 ψψ state must also be an SU (N ) singlet. The different QCD representations for ψψ 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 ψψ has to be an SU (N ) singlet, the new potential will be of the form where V SU (3) is given in Eq. (3.5). The ψψ → 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 where V SU (N ) is given in Eq. (5.2). Note here that ψψ now is an SU (3) singlet. Finally, the ψψ → 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 ψψ 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 The Sommerfeld-corrected annihilation cross sections for all these processes in the s-wave can then be expressed as α N β , (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 α s (μ) = α N (μ) = 0.1. For typical freeze-out velocities v ∼ 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 ψψ → 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].

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 higherorder 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 Sommerfeldcorrection 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 com-bination 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 spinscalar, 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.

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 [44]. In order for the amplitudes to be of the form of Eq. (2.6), we work in the so-called final frame where the momenta of the final state are along the zaxis and the momenta of the initial states are characterized by the angles θ and φ. In order for the polarization vectors of the gluons and of the vector to be well-defined, we compute the amplitudes in the helicity basis based on [54]. In this basis, for the annihilation processes we consider in the main body of the paper the φ dependence of the amplitude is well-known 4 For these processes, we therefore compute the amplitudes for φ = 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 's into quark and gluon pairs are with E = p 2 + m 2 . The spinors for a particle of helicity ± 1 2 and mass m moving in the direction (θ, φ) are The spinors for a particle moving in the opposite direction are obtained in [54] as well. They are very similar to the expressions in Eq. (A.3) and are given in the attached Mathematica notebook [44]. 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 of momentum p 1 are while the longitudinal polarization vector corresponding to a vector of momentum p 1 is Similarly, the transverse polarization vectors corresponding to a final-state gluon of momentum p 4 or to a vector of momentum p 2 are (2) while the longitudinal polarization vector corresponding to a vector of momentum p 2 is The full φ-dependence for polarization vectors can be found in the attached Mathematica notebook [44]. We use [55] 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 where the D s m,λ are the Wigner D-functions. In this formula, we used λ 3 = m 3 and λ 4 = −m 4 for final states, whose momenta are along the z-axis.

B Color decomposition
In this appendix we describe the decomposition of the group structure of the amplitudes discussed in Sect. 3.2 for QCD and in Sect. 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 Sect. 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 R ⊗ R with R = F (fundamental), S ( N (N +1) 2 -dimensional symmetric) and A (adjoint). The decomposition for these channels is given in Eq. (3.3) for QCD and in Eq. (5.2) 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 [38] and use fundamental indices for all representations. To switch between (T a R ) i j where i, j run from 1 to d 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 in A and drop the color dependence of the remaining part in the amplitudes. This method has been put forward already in [17] 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 [56] and ColorMath [57]. 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 Here the indices i, j = 1, . . . , N represent the color of the i, j . For the product of two symmetric representations S the situation is slightly more complicated as one needs to represent each u where u = 1, . . . , 1 2 N (N + 1) with two fundamental indices i, j = 1, . . . , N . We can now write A i j kl = v i j w kl , which now has to be symmetric under the transformations i ↔ j, k ↔ l. Transforming between both representations can be done using the respective Clebsch-Gordan coefficients [58]. The symmetricity representation decomposes as For the adjoint we can write A i j Again we can transform to adjoint indices a = 1, . . . , N 2 −1 by using the respective Clebsch-Gordan coefficients which are obtained directly from the generators of SU (N ). The adjoint decomposes as The above decomposition applies for N ≥ 4, however, in the case of QCD with N = 3 the representation B S does not appear and the symmetric adjoint representation for SU (3) is given by kl . This concludes the decomposition of the amplitudes considered in Sects. 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 Eq. (4.3). 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: By virtue of Eq. (3.2) while inserting the quadratic Casimir invariants from Eq. (3.4) we decompose the QCD potential as .

(B.5)
To decompose the total cross section we write the color part of the amplitude in tensor notation as A i jk lmn = v i jk w lmn with full symmetricity in the upper and lower components separately. After doing the calculation we find [1]   We observe here that in contrast to the correction factors for the triplet, sextet and octet in Eq. (3.18), 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 To decompose the total cross section we write the color part of the amplitude in tensor notation as A 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 3 ⊗ 8 → 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. will have momentum-dependent factors α, β and γ , 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 Sect. 3.3 to each of these terms will give the total analytic Sommerfeldcorrected cross section at a given order in the partial wave expansion. For convenience we present here the color decom-positions of each combination of channels: Here, we consider a model involving a scalar dark matter particle that couples to a Standard Model quark and a new vector-like quark (ψ) which is a triplet under color [37]. In this case, the annihilation of ψ into quark pairs can occur through two independent processes, namely ψψ → qq and ψ ψ → q q plus its conjugate. Annihilation into a quarkanti-quark pair can occur through either an s-channel gluon or a t-channel scalar (the dark matter). ψ 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 ψψ → qq. The corresponding potential has been derived in (3.5) and the resulting Sommerfeld correction for the s-channel gluon exchange has been derived in Sect. 3.2 and presented in Eq. (3.18). Including the new tchannel 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 Sect. B.3 and observe that the t-channel amplitude decomposes as color A t-channel 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 ψ ψ → 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 To decompose the total cross section we write the color part of the initial states in tensor notation as A i j = v i w j and find [3] i j = 1 2 (A i j − A ji ), The total color structure for this process can be written as the sum of the t-channel and u-channel contributions, namely which, beyond the s-wave, contrasts with the result derived in [37] for a similar simplified model.