Studies of a thermally averaged p-wave Sommerfeld factor

TL;DR

This paper examines how attractive SU(3) gauge interactions affect the thermal p-wave annihilation rate of heavy particles, revealing significant enhancement factors up to 100 through theoretical and lattice methods.

Contribution

It provides the first detailed analysis of thermally averaged p-wave Sommerfeld factors including bound state effects for SU(3) interactions.

Findings

Bound states significantly enhance annihilation rates.

Enhancement factors can reach up to approximately 100.

Both perturbative and lattice methods agree on the estimates.

Abstract

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.

Equations

$$\dot{n}+3Hn=-\langle\sigma v\rangle\,\bigl{(}n^{2}-n_{\mbox{\scriptsize eq}}^{2}\bigr{)}\;,$$

$$⟨\sigma v⟩=\frac{{\Gamma }_{\text{mbox}chem}}{2n_{\text{mbox}eq}}.$$

$$⟨\sigma v⟩=4\sum _i{c}_i\frac{⟨{\mathcal{O}}_i⟩}{n_{\text{mbox}eq}^2}.$$

$$\sigma v={\sigma }_{\text{mbox}tree}v\times S(v),$$

$$⟨\sigma v⟩\simeq \frac{{\int }_{\mathbf{v}}\sigma v{e}^{-M_{\text{mbox}\mathrm{kin}}{v}^2/T}}{{\int }_{\mathbf{v}}{e}^{-M_{\text{mbox}\mathrm{kin}}{v}^2/T}}.$$

$${\bar{S}}_i\equiv \frac{⟨{\mathcal{O}}_i⟩/⟨{\mathcal{O}}_i{⟩}_{\text{mbox}tree}}{n_{\text{mbox}eq}^2/(n_{\text{mbox}eq}^2{)}_{\text{mbox}tree}}.$$

$$⟨\sigma v⟩=4\sum _i{c}_i{\bar{S}}_i\frac{⟨{\mathcal{O}}_i{⟩}_{\text{mbox}tree}}{(n_{\text{mbox}eq}^2{)}_{\text{mbox}tree}},$$

$${\mathcal{O}}_p$$

$$\bigl{\langle}{\mathcal{O}}_{p}\bigr{\rangle}_{\mbox{\scriptsize tree}}\;=\;N\int_{{\bf p,q}}\frac{({\bf p}-{\bf q})^{2}}{4M_{\mbox{\scriptsize kin}}^{2}}e^{-(E_{p}+E_{q})/T}\;=\;N\times\frac{3T}{2M_{\mbox{\scriptsize kin}}}\times\biggl{(}\frac{M_{\mbox{\scriptsize kin}}T}{2\pi}\biggr{)}^{3}e^{-{2M_{\mbox{\scriptsize rest}}}/{T}}\;,$$

$$\bigl{(}n_{\mbox{\scriptsize eq}}\bigr{)}_{\mbox{\scriptsize tree}}\;\equiv\;2N\int_{{\bf p}}e^{-E_{p}/T}=2N\,\biggl{(}\frac{M_{\mbox{\scriptsize kin}}T}{2\pi}\biggr{)}^{{\frac{3}{2}\,}}e^{-{M_{\mbox{\scriptsize rest}}}/{T}}\;,$$

$$\frac{⟨{\mathcal{O}}_p{⟩}_{\text{mbox}tree}}{(n_{\text{mbox}eq}^2{)}_{\text{mbox}tree}}=\frac{3T}{8N{M}_{\text{mbox}kin}}.$$

$$E_p+E_q=2M_{\text{mbox}rest}+\frac{k^2}{4M_{\text{mbox}kin}}+E^{\prime },\mathbf{k}\equiv \mathbf{p}+\mathbf{q}.$$

$$\displaystyle\biggl{\{}-\frac{\nabla_{\bf r}^{2}}{M_{\mbox{\scriptsize kin}}}+V_{\mbox{\tiny\rm{$T$}}}(r)-i\Gamma_{\mbox{\tiny\rm{$T$}}}(r)-E^{\prime}\biggr{\}}\,{\bf G}(E^{\prime};{\bf r},{\bf r}^{\prime})$$

$$\displaystyle\lim_{{\bf r},{\bf r}^{\prime}\to{\bf 0}}\mathop{\mbox{Im}}\bigl{[}\nabla_{{\bf r}}\cdot{\bf G}(E^{\prime};{\bf r},{\bf r}^{\prime})\bigr{]}$$

$$⟨{\mathcal{O}}_p⟩$$

$${\bar{S}}_p$$

$$\rho_{\mbox{\scriptsize$p$,tree}}(E^{\prime})\;\equiv\;\frac{N}{M_{\mbox{\scriptsize kin}}^{2}}\int_{{\bf p}}{p}^{2}\pi\,\delta\Bigl{(}E^{\prime}-\frac{p^{2}}{M_{\mbox{\scriptsize kin}}}\Bigr{)}\;=\;\frac{NM_{\mbox{\scriptsize kin}}^{1/2}\theta(E^{\prime})(E^{\prime})^{3/2}}{4\pi}\;.$$

$$\rho_{p}(E^{\prime})=\rho_{\mbox{\scriptsize$p$,tree}}(E^{\prime})\,S_{p}(v)\;,$$

$$S_{p}(v)=S_{s}(v)\,\biggl{(}1+\frac{\alpha^{2}}{4v^{2}}\biggr{)}\;,\quad S_{s}(v)\;\equiv\;\frac{\pi\alpha/v}{1-e^{-\pi\alpha/v}}\;.$$

$$\frac{\rho_{s}(E^{\prime})}{M_{\mbox{\scriptsize kin}}^{2}}\;=\;\frac{N\alpha}{4\pi}\int_{0}^{\infty}\!{\rm d}\rho\,\mathop{\mbox{Im}}\biggl{(}\frac{1}{u_{0}^{2}}\biggr{)}\;,\quad\frac{\rho_{p}(E^{\prime})}{M_{\mbox{\scriptsize kin}}^{2}}\;=\;\frac{N\alpha^{3}}{16\pi}\int_{0}^{\infty}\!{\rm d}\rho\,\mathop{\mbox{Im}}\biggl{(}\frac{1}{u_{0}^{2}}+\frac{36}{u_{1}^{2}}\biggr{)}\;,$$

$$\biggl{\{}\frac{\partial^{2}}{\partial\rho^{2}}-\frac{\ell(\ell+1)}{\rho^{2}}+\frac{E^{\prime}+i\Gamma_{\mbox{\tiny\rm{$T$}}}(r)-V_{\mbox{\tiny\rm{$T$}}}(r)}{M_{\mbox{\scriptsize kin}}\alpha^{2}}\biggr{\}}\,u_{\ell}(\rho)=0\;,$$

$$\int_{\rho_{0}}^{\infty}\!{\rm d}\rho\,\mathop{\mbox{Im}}\biggl{(}\frac{1}{u_{1}^{2}}\biggr{)}=\mathop{\mbox{Im}}\Biggl{\{}\frac{1}{u_{1}(\rho_{0})\bigl{[}u^{\prime}_{1}(\rho_{0})+u_{1}(\rho_{0})\frac{1-ik^{3}\rho_{0}^{3}}{\rho_{0}(1+k^{2}\rho_{0}^{2})}\bigr{]}}\Biggr{\}}\;.$$

$$V_{\mbox{\tiny\rm{$T$}}}(r)=-\frac{g_{s}^{2}C_{\mbox{\tiny\rm{F}}}\exp(-m_{\mbox{\tiny\rm{D}}}r)}{4\pi r}\;,\quad\Gamma_{\mbox{\tiny\rm{$T$}}}(r)=\frac{g_{s}^{2}C_{\mbox{\tiny\rm{F}}}T}{2\pi}\int_{0}^{\infty}\!\frac{{\rm d}z\,z}{(z^{2}+1)^{2}}\biggl{[}1-\frac{\sin(zm_{\mbox{\tiny\rm{D}}}r)}{zm_{\mbox{\tiny\rm{D}}}r}\biggr{]}\;,$$

$${\bar{S}}_p\equiv \frac{P_p/P_p^{\text{mbox}cold}}{(P_1/P_1^{\text{mbox}cold}{)}^2},$$

$$G_{\alpha \gamma ;kl}^{\theta }({\tau }_2,\mathbf{x};{\tau }_1,\mathbf{y})$$

$$G_{\alpha \gamma ;kl}^{\chi }({\tau }_2,\mathbf{x};{\tau }_1,\mathbf{y})$$

$$G^{\chi}(\tau_{2},{\bf x};\tau_{1},{\bf y})=-\bigl{[}G^{\theta}(\tau_{1},{\bf y};\tau_{2},{\bf x})\bigr{]}^{\dagger}\;.$$

$$D_{i}G^{\theta}_{\alpha\gamma;kl}\;\equiv\;\bigl{\langle}(D_{i}\theta)_{\alpha k}(\beta,{\bf x})\,\theta^{\dagger}_{\gamma l}(0,{\bf x})\bigr{\rangle}\;,\quad G^{\theta}_{\alpha\gamma;kl;i}\;\equiv\;\bigl{\langle}\theta_{\alpha k}(\beta,{\bf x})\,(D_{i}\theta)^{\dagger}_{\gamma l}(0,{\bf x})\bigr{\rangle}\;.$$

$$\bigl{(}n_{\mbox{\scriptsize eq}}\bigr{)}_{\mbox{\scriptsize latt}}\;=\;2\mathop{\mbox{Re}}{\rm Tr\,}\bigl{\langle}G^{\theta}(\beta,{\bf 0};0,{\bf 0})\bigr{\rangle}\;.$$

$$P_1$$

$$\langle{\mathcal{O}}_{p}\rangle=\frac{1}{2M_{\mbox{\scriptsize kin}}^{2}}\sum_{i=1}^{3}\mathop{\mbox{Re}}{\rm Tr\,}\bigl{\langle}D_{i}G^{\theta}_{;i}G^{\theta\dagger}\,-\,D_{i}G^{\theta}G^{\theta\dagger}_{;i}\bigr{\rangle}\;.$$

$$P_p$$

Full text

June 2019

{centering}

**Studies of a thermally averaged $p$-wave Sommerfeld factor **

Seyong Kima,b and M. Laineb

${}^{\mbox{\scriptsize a}}$*Department of Physics, Sejong University, Seoul 143-747, South Korea

${}^{\mbox{\scriptsize b}}$*AEC, Institute for Theoretical Physics, University of Bern,

Sidlerstrasse 5, CH-3012 Bern, Switzerland

Abstract

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.

1 Introduction

Inelastic processes between a dilute ensemble of heavy particles moving slowly in a thermal environment are encountered in many physical situations. A classic example is given by nuclear reactions taking place within the electromagnetic plasma of stars [1]. In particle physics, we may consider heavy dark matter particles pair-annihilating into Standard Model particles in the early universe, or a heavy quark and anti-quark pair-annihilating into light quarks and gluons in a quark-gluon plasma generated in heavy ion collision experiments.

The theoretical treatment of slow annihilation processes is facilitated by noting that the average kinetic energy of the annihilating particles is small compared with their rest mass, $Mv^{2}\sim T\ll M$. Such a scale separation permits for a factorized description of annihilation processes in terms of a series of long-distance matrix elements times short-distance Wilson coefficients [2]. In particular, the thermal average of an annihilation rate can be expanded as $\langle\sigma v\rangle=a+b\,\langle v^{2}\rangle+\cdots\;$, where $v$ denotes the relative velocity. The term $a$ is said to originate from “$s$-wave” matrix elements, whereas $b$ may be associated with “$p$-wave” ones.

In the presence of long-range interactions, the coefficients $a$ and $b$ may get large corrections compared with a tree-level treatment. For scattering states, this is known as the “Sommerfeld (-Gamow-Sakharov) effect” [3, 4, 5, 6]. Sommerfeld factors are nowadays routinely included in Boltzmann equations for dark matter pair annihilation (cf., e.g., refs. [7, 8, 9, 10, 11]).

Long-range interactions may also lead to the appearance of bound states in the dark sector, which opens up a fast pair-annihilation channel (cf., e.g., refs. [12, 13]). Bound states are particularly important if the dark sector contains particles charged under QCD, as is the case for instance in a prototypical model in which dark matter is a singlet Majorana fermion and the mediator is a slightly heavier strongly coupled scalar (cf. refs. [14, 15] for reviews).

Recently, we have developed a framework which permits to estimate the thermally averaged pair annihilation rate, including bound-state effects, beyond perturbation theory [16]. The framework can be applied to a number of cosmological models [17], particularly the prototypical framework mentioned above [18, 19], where bound-state effects have been seen to be important from other considerations as well [20, 21, 22, 23].

The purpose of the present work is to extend ref. [16] from the $s$-wave to the $p$-wave case. Even if the $p$-wave contribution is suppressed by $\langle v^{2}\rangle$, its “standard” Sommerfeld enhancement is larger than for $s$-wave [24, 25]. If the coefficients of the $s$-wave operators happen to vanish at leading order, $p$-wave may be the dominant channel [26]. $p$-wave annihilation has also been discussed from astrophysical motivations (cf., e.g., refs. [27, 28, 29]).

This presentation is organized as follows. After outlining the basic setup (cf. sec. 2), we review thermally averaged pair annihilation rates within resummed perturbation theory (cf.  sec. 3). Having introduced the lattice framework (cf. sec. 4), we present and discuss numerical results (cf. sec. 5), and conclude then with a brief summary (cf. sec. 6).

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

[TABLE]

where $H$ is the Hubble rate, $\sigma$ is an annihilation cross section, $v$ is a relative velocity, and $\langle...\rangle$ indicates a thermal average over the momenta of the annihilating particles.

If the dark sector experiences strong interactions, the thermal average $\langle\sigma v\rangle$ may receive large radiative corrections. In order to address these beyond perturbation theory, it was noted in ref. [33] that by linearizing eq. (2.1) close to equilibrium, we may interpret the averaged cross section as being related to a chemical equilibration rate ($\equiv\Gamma_{\mbox{\scriptsize chem}}$),

[TABLE]

Subsequently we can make use of linear response theory in order to relate $\Gamma_{\mbox{\scriptsize chem}}$ to an equilibrium correlator. Furthermore, if we find ourselves in the non-relativistic regime, i.e. with dark matter masses $M\gg\pi 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 [16],

[TABLE]

Here the Wilson coefficients $c_{i}$ and the operators ${\mathcal{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 $\langle...\rangle$.

As is usual for effective field theories, the operators ${\mathcal{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 $\langle\nabla^{2}\rangle/M^{2}\sim\pi T/M\ll 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 [24, 25] and also display bound states, and they thus merit a detailed look.111It has been suggested that, apart from influencing the value of $\langle\sigma v\rangle$, bound states also lead to a modification of the functional form of the part $n^{2}-n_{\mbox{\scriptsize eq}}^{2}$ in eq. (2.1) at late times when $n-n_{\mbox{\scriptsize eq}}\gg n_{\mbox{\scriptsize eq}}$ so that we leave the linear response regime [34]. Furthermore, when $\pi T\ll\Delta E$, where $\Delta 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

[TABLE]

after which thermal averaging is often implemented as

[TABLE]

In reality, vacuum and thermal effects cannot be factorized in this way. Indeed thermal corrections can also modify masses like $M_{\mbox{\scriptsize rest}}$ and $M_{\mbox{\scriptsize 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.

[TABLE]

Here we define $\langle{\mathcal{O}}_{i}\rangle_{\mbox{\scriptsize tree}}$ and $(n_{\mbox{\scriptsize eq}}^{2})_{\mbox{\scriptsize 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}_{\mbox{\scriptsize eq}}$ and $\langle{\mathcal{O}}_{i}\rangle$ by a large amount [1]. As a consequence of this definition, eq. (2.3) can now be re-expressed as

[TABLE]

where the tree-level ratio ${\langle{\mathcal{O}}_{i}\rangle_{\mbox{\scriptsize tree}}}/{(n_{\mbox{\scriptsize eq}}^{2})_{\mbox{\scriptsize 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 [16], these fields are taken to be a spin-${\frac{1}{2}\,}$ particle and antiparticle (that is, heavy quark and antiquark), each with $N\equiv 2N_{\rm c}$ degrees of freedom. However, spin-dependent effects are highly suppressed, so we believe our results to be valid also for spin-[math] particles, such as stops, with the replacement $N\to N_{\rm c}$. The particle and antiparticle fields are denoted by $\theta$ and $\chi$, respectively, and the annihilation operator considered is defined in eq. (3.1).

3 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

[TABLE]

Here $\theta$ and $\chi^{\dagger}$ are annihilation operators for particles and antiparticles, respectively. As the annihilation operators appear on the right, the vacuum state does not contribute to $\langle{\mathcal{O}}_{p}\rangle$.

It is straightforward to evaluate the thermal expectation value of eq. (3.1) in tree-level perturbation theory. We obtain

[TABLE]

where $E_{p}\equiv M_{\mbox{\scriptsize rest}}+p^{2}/(2M_{\mbox{\scriptsize kin}})$ is a non-relativistic energy.222At $T>0$, $M_{\mbox{\tiny\rm{rest}}}$ and $M_{\mbox{\tiny\rm{kin}}}$ do not coincide because of the so-called Salpeter correction to $M_{\mbox{\tiny\rm{rest}}}$, cf., e.g., refs. [1, 35]. Even in vacuum the two can differ if UV regularization does not respect Lorentz invariance, as is the case for instance within the lattice NRQCD setup.

Similarly,

[TABLE]

and correspondingly

[TABLE]

This displays a characteristic $p$-wave velocity suppression by $T/M_{\mbox{\scriptsize kin}}\ll 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

[TABLE]

Moreover, it is useful to resolve $\langle{\mathcal{O}}_{p}\rangle$ into a spectral representation, so that contributions from soft energy scales can be inspected more carefully. A thermal potential V_{\mbox{\tiny\rm{T}}}(r) (cf. eq. (3.16)) is assumed normalized so that \lim_{r\to\infty}V_{\mbox{\tiny\rm{T}}}(r)=0, i.e. $r$-independent thermal corrections, known as the Salpeter correction, have been included in the definition of $M_{\mbox{\scriptsize rest}}$. A vector-like Green’s function is solved for from

[TABLE]

where $\rho_{p}(E^{\prime})$ is a spectral function. Following refs. [16, 17] and carrying out the integral over the center-of-mass momentum ${\bf k}$ (cf. eq. (3.5)), we then get

[TABLE]

Here $\alpha^{2}M_{\mbox{\scriptsize kin}}\ll\Lambda\ll M_{\mbox{\scriptsize rest}}$ is a cutoff restricting the average to the non-relativistic regime. As our masses $M_{\mbox{\scriptsize rest}}$, $M_{\mbox{\scriptsize kin}}$ already include thermal corrections, $n_{\mbox{\scriptsize eq}}=(n_{\mbox{\scriptsize eq}})_{\mbox{\scriptsize 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_{\mbox{\tiny\rm{T}}}(r)\to 0 and \Gamma_{\mbox{\tiny\rm{T}}}(r)\to 0^{+}, eq. (3.6) is easily solved in momentum space, yielding ($p\equiv|{\bf p}|$)

[TABLE]

Inserting into eq. (3.9) and carrying out the integral over $E^{\prime}$ indeed gives unity.

Another limit in which $\rho_{p}(E^{\prime})$ can be determined is a Coulombic potential, namely V_{\mbox{\tiny\rm{T}}}(r)\to-\alpha/r and \Gamma_{\mbox{\tiny\rm{T}}}(r)\to 0^{+}. Parametrizing $E^{\prime}=M_{\mbox{\scriptsize kin}}v^{2}$, the above-threshold solution reads

[TABLE]

where $S_{p}$ is a vacuum $p$-wave Sommerfeld factor, given by (cf., e.g., refs. [24, 25])

[TABLE]

A large enhancement is observed for $v\ll\alpha$, in particular $\lim_{E^{\prime}\to 0}\rho_{p}(E^{\prime})=N\alpha^{3}M_{\mbox{\scriptsize kin}}^{2}/16$. This enhancement originates from an overlap with an $s$-wave radial function (this is explained in footnote 3), and gives the dominant above-threshold contribution to $\bar{S}_{p}$ if T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha^{2}M_{\mbox{\scriptsize 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. [36], and an implementation for the $p$-wave was worked out in ref. [37]. The solutions can be written as333Let us elaborate on the origin of the two parts in $\rho_{p}$. In terms of eigenstates of the operator in eq. (3.6), the $p$-wave solution contains $\nabla\psi({\bf 0})$. In spherical coordinates, writing $\psi=R_{nl}(r)Y_{lm}(\Omega)$, we thus need $R_{nl}^{\prime}(0)$. In a Coulomb potential, $R_{n0}$ has a linear term at small $r$, $R_{n0}(r)=c_{0}+c_{1}r+...$, which leads to an $s$-wave contribution to $\nabla\psi({\bf 0})$, denoted in eq. (3.13) by $u_{0}\sim rR_{n0}$. This is responsible for the dominant term $\sim\alpha^{2}/(4v^{2})$ in eq. (3.12). The second term in $\rho_{p}$ of eq. (3.13) is the “genuine” $p$-wave contribution, originating from the short-distance asymptotics $R_{n1}(r)=d_{0}r+...$.

[TABLE]

where $\alpha\equiv g_{s}^{2}C_{\mbox{\tiny\rm{F}}}/(4\pi)$, $\rho\equiv r\alpha M_{\mbox{\scriptsize kin}}$, and $u_{\ell}$ is a regular solution of the homogeneous equation

[TABLE]

assumed normalized as $u_{\ell}=\rho^{\ell+1}+...$ at short distances. It is sufficient to solve the equation up to some finite $\rho\gg 1$ and attach this to the known asymptotics.444For $\rho\gg 1$, $V_{T}$ vanishes and \Gamma_{\mbox{\tiny\rm{T}}} goes over to a constant, whereby the equation satisfied by the $p$-wave wave function reads $(\partial_{\rho}^{2}-\frac{2}{\rho^{2}}+\hat{E}^{\prime}+i\hat{\Gamma})u_{1}=0$, where $\hat{E}^{\prime}\equiv E^{\prime}/(M_{\mbox{\tiny\rm{kin}}}\alpha^{2})$ and \hat{\Gamma}\equiv\Gamma_{\mbox{\tiny\rm{T}}}(\infty)/(M_{\mbox{\tiny\rm{kin}}}\alpha^{2}). Let us denote $k\equiv\sqrt{\hat{E}^{\prime}+i\hat{\Gamma}}$. Then the general solution reads $u_{1}=C[\sin(k\rho+\delta)/(k\rho)-\cos(k\rho+\delta)]$, where $C,\delta\in\mathbbm{C}$. The function $1/u_{1}^{2}$ is integrable, and subsequently $C\sin(k\rho+\delta)$ and $C\cos(k\rho+\delta)$ can be expressed in terms of $u_{1}(\rho)$ and $u^{\prime}_{1}(\rho)$. We thus obtain

\int_{\rho_{0}}^{\infty}\!{\rm d}\rho\,\mathop{\mbox{Im}}\biggl{(}\frac{1}{u_{1}^{2}}\biggr{)}=\mathop{\mbox{Im}}\Biggl{\{}\frac{1}{u_{1}(\rho_{0})\bigl{[}u^{\prime}_{1}(\rho_{0})+u_{1}(\rho_{0})\frac{1-ik^{3}\rho_{0}^{3}}{\rho_{0}(1+k^{2}\rho_{0}^{2})}\bigr{]}}\Biggr{\}}\;.

(3.15) Setting $\rho_{0}\to\epsilon\equiv 0^{+}$ and recalling $u_{1}(\epsilon)\approx\epsilon^{2}$ yields $\mathop{\mbox{Re}}(k^{3}/9)$, which reproduces eq. (3.10) from eq. (3.13). The part \int_{\rho_{0}}^{\infty}\!{\rm d}\rho\,\mathop{\mbox{Im}}\bigl{(}{1}/{u_{0}^{2}}\bigr{)} of eq. (3.13) yields $\mathop{\mbox{Re}}(k)$ when $\rho_{0}\to\epsilon$, which amounts to $\sim\alpha^{2}\theta(E^{\prime})(E^{\prime})^{1/2}$.

As far as the potential goes, at large separations we make use of a Hard Thermal Loop resummed thermal expression which includes the effects of Debye screening and Landau damping [38, 39, 40],

[TABLE]

where $m_{\mbox{\tiny\rm{D}}}\sim g_{s}T$ is a Debye mass. For numerical estimates we insert 2-loop values of $m_{\mbox{\tiny\rm{D}}}$ and $g_{s}^{2}$ from ref. [41] (the 3-loop level has been reached only for $m_{\mbox{\tiny\rm{D}}}$ [42]). The real part of the potential is interpolated into a vacuum expression at short separations [43, 44], as discussed in ref. [45]. In order to account for the proper kinematics of real processes in a regime beyond which the derivation is formally valid, we also follow the arguments presented in ref. [46] and multiply the imaginary part of the potential by the Boltzmann factor $e^{-|E^{\prime}|/T}$ for $E^{\prime}<0$. Corresponding numerical solutions of the spectral functions $\rho_{s}$ and $\rho_{p}$ are shown in fig. 1.

4 Lattice framework

On the lattice the double ratio in eq. (2.6) is replaced through

[TABLE]

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 that $\bar{S}_{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_{\mbox{\scriptsize 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^{\mbox{\scriptsize th}}$ direction with origin at ${\bf 0}$, by ${\bf i}\equiv a_{s}{\bf e}_{i}$ a displacement in the $i^{\mbox{\scriptsize th}}$ direction by a lattice spacing $a_{s}$, and by $G^{\theta}$, $G^{\chi}$ the propagators

[TABLE]

where $\alpha,\gamma\in\{1,...,N_{\rm c}\}$ are colour indices and $k,l\in\{1,2\}$ are spin indices. Given that $\chi$ represents an antiparticle to $\theta$, the two propagators are related by

[TABLE]

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 $\beta\equiv 1/T$. For taking derivatives of a propagator with respect to the position of a sink or source we introduce a shorthand notation,

[TABLE]

With these propagators, the lattice analogue of $n_{\mbox{\scriptsize eq}}$ reads [16]

[TABLE]

Given that overall normalization cancels out in eq. (4.1), we in practice define $P_{1}$ by dividing $(n_{\mbox{\scriptsize eq}})_{\mbox{\scriptsize latt}}$ by the number of degrees of freedom, viz. $2N$, i.e.

[TABLE]

For the operator in eq. (3.1), Wick contractions yield

[TABLE]

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

[TABLE]

The diagrams illustrate the topology of the contractions.

The lattice framework and the gauge ensemble are as in ref. [16]. The light sector consists of SU(3) gauge theory and $N_{\rm f}=2+1$ flavours of vectorlike fermions transforming in the fundamental representation. The parameters of the action were tuned in refs. [47, 48]. Denoting by $\Lambda$ a scale parameter (cf. ref. [49] for a review), the lightest pseudoscalar mesons have masses $1.2\Lambda$ and $1.5\Lambda$, respectively, the latter for the mesons involving one quark of the third flavour. The lattice is anisotropic, with $a_{s}/a_{\tau}\approx 3.5$, where the spatial lattice spacing is $a_{s}\approx 0.21\Lambda^{-1}$. The spatial extent of the box is $L=24a_{s}$. The system is put at a finite temperature by tuning $N_{\tau}$, i.e. the number of temporal lattice sites, so that $T=\frac{1}{N_{\tau}a_{\tau}}$. The system has a (pseudo)critical temperature at $T_{\rm c}\approx 0.54\Lambda$ [50]. Thermal properties of the system were studied in ref. [51]. We vary $T=(0.95...1.9)T_{\rm c}$ and, setting $M_{\mbox{\scriptsize kin}}=14\Lambda$, can hence access values $M_{\mbox{\scriptsize kin}}/T\sim 14...28$, a reasonable range in view of dark matter freeze-out computations.

5 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 (the errors shown for the lattice results are statistical only). For the $p$-wave, shown in fig. 2(right), we find surprisingly good qualitative agreement, indicating an enhancement factor $\sim 100$ at the lowest temperature. We note that the system is in a confined phase for M/T\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}26.

For the $s$-wave, shown in fig. 2(left), the discrepancy between the perturbative and lattice results is rather substantial.555In ref. [16], the perturbative values were noticeably larger, and the agreement looked better. There are two reasons for this: in ref. [16] 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.

In fact, naively $\bar{S}_{p}>\bar{S}_{s}$ (cf. eq. (3.12)), whereas on the lattice $\bar{S}_{s}$ clearly exceeds $\bar{S}_{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 ($|\psi({\bf 0})|^{2}$, $|\nabla\psi({\bf 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. [52]).

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 $\alpha_{s}\simeq 0.3...0.6$ for $M_{\mbox{\scriptsize kin}}/T\simeq 10...30$. The grey bands in fig. 2 originate from the variation of a thermal $\alpha_{s}$ [41] as we change the renormalization scale within a factor ${\frac{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 of $\bar{S}_{s}$ is influenced by above-threshold scattering states, i.e. tree-level processes, which are insensitive to $\alpha_{s}$. If we artificially increase $\alpha_{s}$ by a factor two, into the range $0.6...1.2$, then $\bar{S}_{s}$ increases by a factor $3...20$, improving the agreement, however $\bar{S}_{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, including the need for a careful power counting concerning which resummations are necessary in the various temperature and mass ranges of interest.

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.

6 Conclusions

Building upon a framework developed in ref. [16], 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 $\sim 100$ are possible (cf. fig. 2).

Within naive perturbation theory, $\bar{S}_{p}>\bar{S}_{s}$ (cf. eq. (3.12)), but on the lattice we find $\bar{S}_{s}>\bar{S}_{p}$ (cf. fig. 2). We may speculate that the large non-perturbative increase of $\bar{S}_{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\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}1 TeV, we normally find ourselves in the regime $T\gg\Lambda$, which implies that $\alpha$ is smaller than in our study. However, as indicated by eq. (3.12), the magnitude of the averaged Sommerfeld factors depends on the ratio $\sim\langle\alpha/v\rangle\sim\sqrt{\alpha^{2}M_{\rm kin}/T}$. Therefore large averaged Sommerfeld factors are found at least in the regime $M_{\rm kin}/T\gg 100$, relevant for late-time pair annihilations. Because of the smaller $\alpha$, higher order corrections should be smaller than in our study. The fact that we find qualitative resemblances even in our fig. 2, then suggests that resummed perturbative estimates should be conservative in that case. For $M_{\mbox{\scriptsize kin}}/T\in(10,1000)$, resummed perturbative values of $\bar{S}_{s}$ from ref. [18] can be found on the web site http://www.laine.itp.unibe.ch/sommerfeld, and we have now added corresponding results for $\bar{S}_{p}$ there.

Acknowledgements

We thank the FASTSUM collaboration for providing the gauge configurations used in our measurements. M.L. thanks Maxim Khlopov for interesting correspondence. S.K. was supported by the National Research Foundation of Korea under grant No. 2018R1A2A2A05018231 funded by the Korean government (MEST) and in part by NRF-2008-000458. M.L. was supported by the Swiss National Science Foundation (SNF) under grant 200020-168988.