Studies of a thermally averaged p-wave Sommerfeld factor

Thermal pair annihilation of heavy particles, such as dark matter or its co-annihilation partners, can be strongly influenced by attractive interactions. We investigate the case that pair annihilation proceeds through a velocity-suppressed $p$-wave operator, in the presence of an SU(3) gauge force. Making use of a non-relativistic effective theory, the thermal average of the pair-annihilation rate is estimated both through a resummed perturbative computation and through lattice simulation, in the range $M/T \sim 10 ... 30$. Bound states contribute to the annihilation process and enhancement factors of up to $\sim 100$ can be found.


Basic setup
Denoting by n the dark matter number density, and assuming that there is a discrete quantum number which prohibits dark matter from decaying, its cosmological evolution is normally described by the Lee-Weinberg equation [28][29][30], n + 3Hn = − σv n 2 − n 2 eq , (2.1) where H is the Hubble rate, σ is an annihilation cross section, v is a relative velocity, and ... indicates a thermal average over the momenta of the annihilating particles.
If the dark sector experiences strong interactions, the thermal average σv may receive large radiative corrections. In order to address these beyond perturbation theory, it was noted in ref. [31] that by linearizing eq. (2.1) close to equilibrium, we may interpret the averaged cross section as being related to a chemical equilibration rate (≡ Γ chem ), Subsequently we can make use of linear response theory in order to relate Γ chem to an equilibrium correlator. Furthermore, if we find ourselves in the non-relativistic regime, i.e. with dark matter masses M ≫ πT , then the annihilations can be described by local operators [1], similar to those found in the NRQCD context [2]. Then the equilibrium correlators can be reduced to thermal expectation values of the annihilation operators [14], Here the Wilson coefficients c i and the operators O i can be taken over from a vacuum computation, capturing the contribution of "hard scales" to the annihilation process, whereas the influence of the "soft scales" resides within the thermal expectation value ... . As is usual for effective field theories, the operators O i can be organized as an expansion in 1/M 2 . The leading terms, called s-wave operators, do not contain derivatives and are suppressed by 1/M 2 . At the next order, operators appear which contain two spatial derivatives and which are correspondingly suppressed by 1/M 4 . Given that ∇ 2 /M 2 ∼ πT /M ≪ 1, the p-wave operators are normally strongly suppressed compared with the s-wave operators. However, p-wave operators may experience relatively speaking larger enhancements from interactions [22,23] and also display bound states, and they thus merit a detailed look. 1 1 It has been suggested that, apart from influencing the value of σv , bound states also lead to a modification of the functional form of the part n 2 − n 2 eq in eq. (2.1) at late times when n − n eq ≫ n eq so that we leave the linear response regime [32]. Furthermore, when πT ≪ ∆E, where ∆E is a binding energy, bound states fall out of chemical equilibrium, and should be added as separate variables in the set of rate equations.
The way that interactions modify the annihilation process can be parametrized through "Sommerfeld factors". In vacuum, the Sommerfeld factor for an annihilation from unbound states is defined by writing after which thermal averaging is often implemented as In reality, vacuum and thermal effects cannot be factorized in this way. Indeed thermal corrections can also modify masses like M rest and M kin , and open up new channels not present in vacuum, like scatterings off light plasma particles. A proper definition of a thermally averaged Sommerfeld factor can be given for the combination appearing in eq. (2.3) and for each operator separately, viz.
Here we define O i tree and (n 2 eq ) tree as tree-level quantities. The rationale of the double ratio in eq. (2.6) is that it removes effects not only from the tree-level scattering process but also from "trivial" corrections to the rest mass, which affect n 2 eq and O i by a large amount [1]. As a consequence of this definition, eq. (2.3) can now be re-expressed as where the tree-level ratio O i tree /(n 2 eq ) tree is dimensionless and has a simple expression, for instance as given in eq. (3.4) for the operator in eq. (3.1).
To be concrete, we consider a theory with heavy particles charged under the fundamental and antifundamental representation of SU(3). Following the original inspiration from QCD [14], these fields are taken to be a spin-1 2 particle and antiparticle (that is, heavy quark and antiquark), each with N ≡ 2N c degrees of freedom. However, spin-dependent effects are highly suppressed, so we believe our results to be valid also for spin-0 particles, such as stops, with the replacement N → N c . The particle and antiparticle fields are denoted by θ and χ, respectively, and the annihilation operator considered is defined in eq. (3.1).

Perturbative considerations
Assuming that the overall scaling of the annihilation operators as 1/M 2 has been incorporated into the coefficients c i in eq. (2.3), the p-wave operator that we consider is defined as Here θ and χ † are annihilation operators for particles and antiparticles, respectively. As the annihilation operators appear on the right, the vacuum state does not contribute to O p . It is straightforward to evaluate the thermal expectation value of eq. (3.1) in tree-level perturbation theory. We obtain and correspondingly This displays a characteristic p-wave velocity suppression by T /M kin ≪ 1.
In order to determine the perturbative value of the averaged Sommerfeld factor of eq. (2.6), it is helpful to go over into center-of-mass coordinates, defined as Moreover, it is useful to resolve O p into a spectral representation, so that contributions from soft energy scales can be inspected more carefully. A thermal potential V T (r) (cf. eq. (3.16)) is assumed normalized so that lim r→∞ V T (r) = 0, i.e. r-independent thermal corrections, known as the Salpeter correction, have been included in the definition of M rest . A vector-like Green's function is solved for from where ρ p (E ′ ) is a spectral function. Following refs. [14,15] and carrying out the integral over the center-of-mass momentum k (cf. eq. (3.5)), we then get Here α 2 M kin ≪ Λ ≪ M rest is a cutoff restricting the average to the non-relativistic regime. As our masses M rest , M kin already include thermal corrections, n eq = (n eq ) tree within our approximation, so that eq. (3.9) is obtained by dividing eq. (3.8) through (3.2). Let us crosscheck that eqs. (3.6)-(3.9) are correct at tree-level. Setting V T (r) → 0 and Γ T (r) → 0 + , eq. (3.6) is easily solved in momentum space, yielding (p ≡ |p|) Inserting into eq. (3.9) and carrying out the integral over E ′ indeed gives unity. Another limit in which ρ p (E ′ ) can be determined is a Coulombic potential, namely V T (r) → −α/r and Γ T (r) → 0 + . Parametrizing E ′ = M kin v 2 , the above-threshold solution reads where S p is a vacuum p-wave Sommerfeld factor, given by (cf., e.g., refs. [22,23]) . This enhancement originates from an overlap with an s-wave radial function, and gives the dominant above-threshold contribution toS p if T < ∼ α 2 M kin .
A general numerical method to find the solution of the s-wave analogues of eqs. (3.6) and (3.7) was presented in ref. [34], and an implementation for the p-wave was worked out in ref. [35]. The solutions can be written as 3 where α ≡ g 2 s C F /(4π) represents an in principle arbitrary choice of a coupling constant, ρ ≡ rαM kin , and u ℓ is a regular solution of the homogeneous equation

Lattice framework
On the lattice the double ratio in eq. (2.6) is replaced through where P 1 and P p are expectation values to be specified presently (cf. eqs. (4.7) and (4.9)). The superscript "cold" indicates a measurement with all link matrices set to unity; this is an implementation of the "tree-level" prescription of perturbation theory. The division by the respective cold measurement implies thatS p deviates from unity only through the effect of gauge interactions. The normalization by P 2 1 furthermore implies that modifications of the rest mass by gauge interactions are cancelled, an effect which is linearly divergent in lattice regularization and strongly influences n eq (cf. eq. (3.3)).
For a lattice measurement, we choose a simple first-order discretization of the covariant derivatives in eq. (3.1). We denote by U i a link in the i th direction with origin at 0, by i ≡ a s e i a displacement in the i th direction by a lattice spacing a s , and by G θ , G χ the propagators where α, γ ∈ {1, ..., N c } are colour indices and k, l ∈ {1, 2} are spin indices. Given that χ represents an antiparticle to θ, the two propagators are related by Because non-relativistic particles move in the positive time direction only, a non-zero contraction may necessitate propagating across the imaginary time interval, whose extent is β ≡ 1/T . For taking derivatives of a propagator with respect to the position of a sink or source we introduce a shorthand notation, With these propagators, the lattice analogue of n eq reads [14] n eq latt = 2 Re Tr G θ (β, 0; 0, 0) . (4.6) Given that overall normalization cancels out in eq. (4.1), we in practice define P 1 by dividing (n eq ) latt by the number of degrees of freedom, viz. 2N , i.e.
For the operator in eq. (3.1), Wick contractions yield Replacing covariant derivatives by discrete lattice derivatives, and choosing again a convenient normalization, whose effects cancel out in eq. (4.1), we are led to define The diagrams illustrate the topology of the contractions.
The lattice framework and the gauge ensemble are as in ref. [14]. The light sector consists of SU(3) gauge theory and N f = 2 + 1 flavours of vectorlike fermions transforming in the fundamental representation. The parameters of the action were tuned in refs. [45,46]. Denoting by Λ a scale parameter [47], the lightest pseudoscalar mesons have masses 1.2Λ and 1.5Λ, respectively, the latter for the mesons involving one quark of the third flavour. The lattice is anisotropic, with a s /a τ ≈ 3.5, where the spatial lattice spacing is a s ≈ 0.21Λ −1 . The spatial extent of the box is L = 24a s . The system is put at a finite temperature by tuning N τ , i.e. the number of temporal lattice sites, so that T = 1 Nτ aτ . The system has a (pseudo)critical temperature at T c ≈ 0.54Λ [48]. Thermal properties of the system were studied in ref. [49]. We vary T = (0.95...1.9)T c and, setting M kin = 14Λ, can hence access values M kin /T ∼ 14...28, a reasonable range in view of dark matter freeze-out computations.

Numerical results and their uncertainties
Perturbative results for thermally averaged Sommerfeld factors from sec. 3 and lattice results from sec. 4 are compared with each other in fig. 2. For the p-wave, shown in fig. 2(right), we find surprisingly good qualitative agreement, indicating an enhancement factor ∼ 100 at the lowest temperature. We note that the system is in a confined phase for M/T > ∼ 26.
For the s-wave, shown in fig. 2(left), the discrepancy between the perturbative and lattice results is rather substantial. 5 In fact, naivelyS p >S s (cf. eq. (3.12)), whereas on the latticeS s clearly exceedsS p . In this context we note that physically, the thermally averaged Sommerfeld factors are sensitive both to energy levels and to the corresponding "overlaps", or wave functions at origin (|ψ(0)| 2 , |∇ψ(0)| 2 ). For another observable in a similar temperature range, it has been found that while for energy levels there is fair agreement, lattice and perturbative overlaps show substantial discrepancies (cf. fig. 6 in ref. [50]). Let us discuss possible reasons for the discrepancy. Starting with the perturbative side, we are quite close to the confined phase and correspondingly our effective coupling is large, varying in the range α s ≃ 0.3...0.6 for M kin /T ≃ 10...30. The grey bands in fig. 2 originate from the variation of a thermal α s [39] as we change the renormalization scale within a factor 1 2 ...2. In the s-wave case, the corresponding error band looks quite narrow. The reason is that in this parameter range the value ofS s is influenced by above-threshold scattering states, i.e. tree-level processes, which are insensitive to α s . If we artificially increase α s by a factor 5 In ref. [14], the perturbative values were noticeably larger, and the agreement looked better. There are two reasons for this: in ref. [14] we used the larger 1-loop thermal coupling, and most importantly the Salpeter correction (thermal shift of the threshold location to smaller energies) was included in the Sommerfeld factor on the perturbative side. The latter has now be excluded from the definition of the Sommerfeld factor through eq. (2.6) on both the perturbative and lattice side, so we believe the comparison to be fairer. two, into the range 0.6...1.2, thenS s increases by a factor 3...20, improving the agreement, howeverS p increases simultaneously by a factor 4...70, spoiling the agreement on that side. In principle a possible way to reduce these uncertainties would be a systematic higher-order computation, however it represents a daunting task because of the resummations needed.
On the lattice side, no infinite-volume or continuum extrapolation was carried out. A box of a finite size influences the spectrum of scattering states, and given that scattering states contribute to the pair annihilation process, this might imply the presence of finite-volume effects. If the system has tightly bound states, whose Bohr radius is not much larger than the lattice spacing, there may also be large discretization effects. In order to check whether the lattice results are plagued by finite-volume or discretization artifacts, additional sets of simulations are needed, requiring a major computational effort beyond our resources.

Conclusions
Building upon a framework developed in ref. [14], we have estimated the thermally averaged p-wave Sommerfeld factor associated with a particular annihilation channel (cf. eq. (3.1)), both through a resummed perturbative (cf. sec. 3) and through a lattice computation (cf. sec. 4). Both methods suggest that large enhancement factors ∼ 100 are possible (cf. fig. 2).
Within naive perturbation theory,S p >S s (cf. eq. (3.12)), but on the lattice we find S s >S p (cf. fig. 2). We may speculate that the large non-perturbative increase ofS s is due to more prominent bound-state effects in the s-wave, however systematic uncertainties may also play a role (cf. sec. 5), an effect which can hopefully be clarified through future work.
In cosmological applications, with M > ∼ 1 TeV, we normally find ourselves in the regime T ≫ Λ. Then the coupling is weaker than in our study, and averaged Sommerfeld factors are smaller, unless M kin /T ≫ 100. For M kin /T ∈ (10, 1000), values ofS s from ref. [16] can be found on the web site http://www.laine.itp.unibe.ch/sommerfeld, and we have now added resummed perturbative values ofS p there. Based on fig. 2, we may expect such estimates to be conservative, so that they can also be used for practical applications.