Rapid thermal co-annihilation through bound states in QCD

The co-annihilation rate of heavy particles close to thermal equilibrium, which plays a role in many classic dark matter scenarios, can be"simulated"in QCD by considering the pair annihilation rate of a heavy quark and antiquark at a temperature of a few hundred MeV. We show that the so-called Sommerfeld factors, parameterizing the rate, can be defined and measured non-perturbatively within the NRQCD framework. Lattice measurements indicate a modest suppression in the octet channel, in reasonable agreement with perturbation theory, and a large enhancement in the singlet channel, much above the perturbative prediction. The additional enhancement is suggested to originate from bound state formation and subsequent decay. Making use of a Green's function based method to incorporate thermal corrections in perturbative co-annihilation rate computations, we show that qualitative agreement with lattice data can be found once thermally broadened bound states are accounted for. We suggest that our formalism may also be applicable to specific dark matter models which have complicated bound state structures.


Introduction
Heavy particle co-annihilation represents a subtle problem, because the slowly moving particles experience strong initial state effects before the final inelastic process. Reversing the time direction, the same can be said about pair creation. Nevertheless, within two-body quantum mechanics in a system with a Coulomb potential, this physics was understood long ago [1,2]. It also plays an essential role for heavy particle pair creation in QCD [3]. The enhancement or suppression of the annihilation rate, depending on whether the Coulomb interaction is attractive or repulsive, is generically referred to as the Sommerfeld effect.
If the co-annihilating particles are part of a thermal medium, the Sommerfeld effect plays a role in a modified form. In this case the velocities of the annihilating particles are not fixed, but come from a statistical distribution. In particular, if we assume the heavy particles to be in kinetic equilibrium, the velocities are distributed according to the Boltzmann weight. Then we speak of thermally averaged Sommerfeld factors. The thermally averaged Sommerfeld factors are relevant for cosmology, particularly for determining the abundance of heavy weakly interacting dark matter particle species [4][5][6][7], which decouple deep in the non-relativistic regime. A recent example of an embedding of the corresponding Sommerfeld factors in a realistic setting can be found in ref. [8].
The starting point of the present paper is the observation that physics formally similar to dark matter co-annihilation takes place with heavy quarks and antiquarks in hot QCD at a temperature of a few hundred MeV. Going beyond a leading-order computation based on Boltzmann equations [9][10][11], the rate of their chemical equilibration can be defined on a general level [12], and can be shown to be sensitive to the physics of the thermal Sommerfeld effect [13]. Even though the associated time scale appears to be too long to play a practical role within the life-time ∼ 10 fm/c of a fireball generated in a heavy-ion collision experiment, this analogy nevertheless means that established methods of QCD, including those of lattice QCD, can be used to investigate the problem. In particular, in ref. [13] a strategy was outlined for implementing in Euclidean spacetime the absorptive parts of 4-quark operators, which can be used for describing the decays of quark-antiquark states within the NRQCD effective theory framework [14,15].
The purpose of the present paper is to realize this proposal. A main part is to develop further the theoretical formulation, showing that the measurement can be reduced to 2point functions and that analytic continuation back to Minkowskian spacetime poses no problem. We also carry out an exploratory lattice study, finding an intriguing pattern with an enhancement in the singlet channel much larger than predicted by the standard formulae used in the literature for incorporating thermal Sommerfeld enhancement. An explanation for this observation in terms of bound-state physics is put forward.
The plan of this paper is the following. After defining the basic observables in sec. 2, we show in sec. 3, through a spectral representation and a canonical analysis, how their measurement can be related to purely static 2-point correlation function ratios. The perturbative evaluation of these correlators is discussed in sec. 4, and corresponding lattice measurements are reported in sec. 5. A concluding discussion and an outlook are offered in sec. 6. Readers not interested in the details of the theoretical and lattice analyses are suggested to consult sec. 2.1 and subsequently proceed to secs. 4 and 6.

Physics background
In order to outline the physics problem in a simple setting, we consider heavy quarks and antiquarks, with a pole mass M , placed in a heat bath at a temperature T ≪ M . We stress that even though we make use of QCD terminology, the discussion until eq. (2.4) applies rather generally, with the gauge group replaced as appropriate. In the QCD context the only specific assumption made is that we consider time scales up to some thousands of fm/c, so that the weak decays of the heavy quarks can be omitted; this can be viewed as the analogue of "R-parity conservation" in supersymmetric theories. Thereby the only number-changing reactions are pair annihilations and creations.
In the heavy-quark limit, QCD has an "emergent" symmetry, in that the quark and antiquark numbers are conserved separately. This conservation is only violated by higherdimensional operators suppressed by 1/M . Omitting such operators it is possible to define separate distribution functions for heavy quarks and antiquarks. We assume that both are close to kinetic and chemical equilibrium, whereby the distribution functions have the forms f p =f p ≈ 2N c exp(−E p /T ), where E p ≡ p 2 + M 2 and 2N c counts the spin and colour degrees of freedom. Here p ≡ |p| is the spatial momentum with respect to the heat bath.
Consider now the total number density, n = p (f p +f p ) + n bound , where n bound denotes the density of quark-antiquark bound states. The equilibrium value of the total density reads n eq ≈ 4N c p exp(−E p /T )+O(e −2M/T ), where the last term indicates the bound state contribution. If the equilibrium conditions are evolving, for instance through a Hubble expansion characterized by a rate H, then the system attempts to adjust its number density to this change. Within a Boltzmann equation approach the evolution equation has the form [16] (∂ t + 3H) n ≈ −c (n 2 − n 2 eq ) . (2.1) If we linearize the right-hand side around equilibrium, this can be rewritten as The coefficient Γ chem = 2 c n eq is called the chemical equilibration rate. It tells how efficiently the system is able to re-adjust its density towards the evolving n eq , and encodes the effects of the microscopic processes which can change the quark and antiquark number densities, notably pair creations and annihilations. Note that, unlike eq. (2.1), the form of eq. (2.2) represents a general linearization and is thus valid beyond the Boltzmann approach. Let us stress that Γ chem describes the slow evolution of a number density, averaged over a large volume. The associated frequency scale is ω ∼ Γ chem ∼ T exp(−M/T ) ≪ T . This slow evolution is caused by infinitely many rare individual processes. The individual processes are associated with pair creations and annihilations and carry a large energy, E ∼ 2M ≫ T . We are not studying these individual reactions separately; in fact the energy scale 2M can be integrated out whereby the system can be described by an effective theory known as the NRQCD [14,15], capturing physics at scales |E − 2M | ≪ M . Now, it is important to realize that even though bound states give a strongly suppressed contribution to n eq , they give an enhanced contribution to Γ chem . Indeed, as already mentioned, in n eq ∼ e −M/T the bound state contribution is n bound ∼ e −2M/T . But Γ chem originates from reactions where a quark and antiquark come together, i.e. |∂ t n| ∼ e −2M/T . In a bound state the quark and antiquark are "already" together, and with a less suppressed Boltzmann weight, because of a binding energy ∆E > 0. Therefore the decay rate from bound states is If T < ∼ ∆E ∼ M α 2 s , this contribution dominates the annihilation rate. This leads to a physical picture in which chemical equilibration proceeds via a "two-stage" process: in thermal equilibrium, only a small fraction of quarks and antiquarks form bound states. But it is those bound states which are most efficiently depleted as the temperature decreases. Bound state formation itself is a fast process, without any exponential suppression factors.
Even though the potential significance of bound states has been mentioned in many cosmological studies, such as refs. [17,18], they have not established themselves as a standard ingredient in quantitative estimates of the thermal freeze-out process. One good reason is certainly that bound states only exist in very specific weakly interacting models. But there is also another issue, namely that the computations are usually based on the Boltzmann equation, eq. (2.1), where c = σv is a thermally averaged annihilation cross section. Thereby, if not sufficient care is taken, the basic assumption that the system can be described by dilute single-particle on-shell states may be inadvertently built into the formalism from the beginning. Recently the importance of bound-state effects has been stressed in, e.g., refs. [19][20][21][22][23], which also suggested various ways to include them. The simplest possibility is to add a bound state phase space distribution as an additional variable in a set of Boltzmann equations (even though it is not clear whether a thermally broadened resonance can be accurately treated as an on-shell degree of freedom). Here we follow a different avenue, aiming eventually to offer for a framework permitting to scrutinize the accuracy of the Boltzmann approach. For this purpose, we take the general linearization of eq. (2.2) as a starting point.

Basic definitions
We next recall how Γ chem can be defined, through a linear response type analysis, as a transport coefficient [12], and how its subsequent evaluation reveals the presence of a Sommerfeld effect [13]. We only give the main steps, referring to refs. [12,13] for details.
Let ψ denote the Dirac spinor of a heavy quark. Making use of a representation in which γ 0 = diag(½ 2×2 , −½ 2×2 ), ψ can be expressed in terms of spinors which have two non-zero components: In the following the vanishing components are omitted. The spinor θ can be associated with a quark state, and the conjugate of the spinor χ with an antiquark state. In particular the energy density carried by heavy quarks and antiquarks, which for large M is equal to their rest mass times their number density, can be described through where θ αi annihilates a quark of colour α and spin i, χ * αi does the same for an antiquark, and M ≡ M rest is the heavy quark rest mass. 1 Note that we treat eq. (2.5), without any terms suppressed by 1/M , as a definition; up to a normalization, it represents a Noether charge density related to the NRQCD Lagrangian.
We are concerned with the imaginary-time 2-point correlator of the operator H, Given that H is a conserved charge density, ∆(τ ) is constant within the standard NRQCD theory [14]. In full QCD, however, the number of heavy quarks is not conserved, because heavy quarks and antiquarks can pair annihilate into gluons and light quarks. These processes can be described by adding 4-quark operators to the NRQCD Lagrangian [15]. The 4-quark operators can be classified according to the gauge and spin quantum numbers of the 2-quark "constituent operators" of which they are composed (the 4-quark operators themselves are gauge singlets and Lorentz scalars). We carry out the main discussion in terms of the "singlet" operator, 1 In perturbative considerations, particularly in sec. 4, we make no difference between rest and kinetic masses, assuming that both correspond to a pole mass. It is worth noting that the precise definition in eq. (2.5) is irrelevant since the mass cancels in eq. (2.12); only the kinetic mass is relevant for the Sommerfeld factors to be defined presently. However the rest mass strongly affects the quark number susceptibility χ f , appearing e.g. in eqs. (2.17) and (2.20), and it also appears trivially in eqs. (2.7) and (2.9). whose coefficient has an "absorptive" imaginary part at O(α 2 s ): For future reference let us also write down one of the octet operators, where T a are generators of SU(N c ), normalized as Tr (T a T b ) = δ ab /2. We work to leading non-trivial order in α s (2M ) whereby O(α 3 s ) corrections to the coefficients (cf. e.g. ref. [24]) as well as operators suppressed by the relative velocity v 2 ∼ α 2 s (cf. e.g. ref. [25]) are omitted. If the singlet operator is added to the NRQCD Lagrangian, and the correlator of eq. (2.6) is computed within this theory, then ∆(τ ) is no longer a constant. We assume that the correlator is computed to first order as an expansion in Im and a corresponding spectral function is determined, , then a "transport coefficient" can be defined as Note that the transport coefficient is extracted from a 1/ω tail of the spectral function, rather than from a linear slope as is sometimes the case; the reason is that the interaction responsible for the process considered has been treated as an insertion. The heavy quark chemical equilibration rate is subsequently obtained from [12] where χ f is the quark-number susceptibility related to the heavy flavour. Note that M appears linearly in H (cf. eq. (2.5)) and quadratically in Ω chem (cf. eq. (2.6)) and therefore drops out in this ratio. The singlet operator in eq. (2.7) mediates decays which experience the so-called Sommerfeld enhancement [1][2][3]. Within perturbation theory, omitting bound state effects, the quarknumber susceptibility χ f and the singlet contribution to Ω chem , denoted by δ 1 Ω chem , read [13] whereS 1 is a thermal Sommerfeld factor, related to a vacuum Sommerfeld factor S 1 bȳ Here E rel ≡ M v 2 , with v = |v|, is the energy of the relative motion of the annihilating particles. The normalization is such that without any resummation,S 1 = S 1 = 1. When higher-order perturbative corrections to S 1 are considered, there is particular subseries of them which proceeds in powers of α s /v. These are summed to all orders into the Sommerfeld factor S 1 . After the thermal average in eq. (2.15), v is parametrically of order T /M . Therefore the thermal Sommerfeld effect is important for T < ∼ α 2 s M , whereas it evaluates to unity for T ≫ α 2 s M . Inspired by eqs. (2.13) and (2.14), we may define a non-perturbative Sommerfeld factor as where δ 1 Ω chem is to be evaluated with the operator in eq. (2.7) inserted to first order. Note that the coefficient Im f 1 ( 1 S 0 ) drops out in this ratio because δ 1 Ω chem is linear in it. With this definition the chemical equilibration rate from eq. (2.12) can be expressed as Thereby the physical result has factorized into a high-energy part Im f 1 ( 1 S 0 ) as well as ingredients which can be addressed non-perturbatively. These areS 1 , which can be extracted from a dimensionless ratio of correlation functions (cf. eq. (3.17)), and the susceptibility χ f , which is a standard observable measuring essentially the equilibrium density n eq . The appearance of 1/N c in eq. (2.17) is a consequence of including N c in χ f . For the octet channel, the analogue of eq. (2.14) reads [13] where the perturbative value ofS 8 is given by an average like in eq. (2.15). We now define the non-perturbative value through For the chemical equilibration rate this leads to (2.20) The colour factors follow from the definition of the octet operator in eq. (2.9) and from the inclusion of one N c in χ f . In QCD there is also another octet operator at O(α 2 s , v 0 ), whose coefficient Im f 8 ( 3 S 1 ) is proportional to the number of light flavours N f [15]; for brevity we concentrate on spin-independent operators here but do not expect any qualitative changes from spin-dependent terms.

Reduction to a static 2-point correlator
When the correlator in eq. (2.6) is evaluated by inserting the operator of eq. (2.7) to first order, we are faced with a 3-point function. We start by recalling that in the imaginary-time formalism the 3-point function splits up into a 4-point function. However, we subsequently show that the 4-point function can be reduced to a simple 2-point function, if we are only interested in extracting the transport coefficient Ω chem and not the full spectral shape.

Absorptive parts of 4-quark operators in imaginary time
The imaginary parts of the 4-fermion vertices in NRQCD represent pair annihilations of heavy quarks and antiquarks. Going to the center-of-mass frame, the vertex is a function of the total energy E of the annihilating pair. Expanding the energy dependence in powers of where v denotes the velocity of q in the qq rest frame, the leading term is a constant, i.e. of O(v 0 ). We restrict to the leading order in the non-relativistic expansion and thereby only consider the constant term.
Even though the term is constant, the corresponding vertex is non-local in imaginary time [12]. In general, the imaginary part corresponds to a spectral function, or cut, of a certain Green's function. More specifically, the 4-fermion vertex can be viewed as a limit of a correlator of two quark-antiquark operators. Given a spectral function for the latter, ρ(E), the corresponding imaginary-time correlator is given by where f B is the Bose distribution. Given that we know ρ only for E ≈ 2M , we cut off the low-E contribution by a mass scale Λ, lying parametrically in the range T ≪ Λ ≪ 2M . Denoting the constant value by ρ 0 , we may then estimate where we approximated |τ | ∼ 1/(2M ) ≪ β, exp(−Λ/T ) ≪ 1, and |τ |Λ ∼ Λ/(2M ) ≪ 1. Therefore, rather than obtaining δ(t) as would be the case in real time, we now get a nonlocality 1/|τ |. Specifically, the imaginary part of the singlet operator of eq. (2.7) can be expressed as a real effective operator in a Euclidean action, We remark that for the argumentation in eq. (3.2) it was essential that 1/|τ | originated from the contribution of large energies and corresponds therefore to a short separation. In eq. (3.3) this restriction has been lifted; this is not a concern because, as our subsequent analysis shows, the contribution still effectively emerges from short separations (cf. eq. (3.14)). A related point is that it is not clear how eq. (3.3) should be defined at the contact point τ 1 = τ 2 . Fortunately, the contact point is naturally regularized by the considerations to follow, so for the moment we simply treat eq. (3.3) as a formal construction.

Canonical analysis and analytic continuation
and making use of translational invariance in order to adjust the spatial coordinates, eqs. (2.6) and (3.3) imply that we are faced with the 4-point correlator We proceed to giving a canonical interpretation to eq. (3.5), showing that the 4-point correlator reduces to a 2-point correlator once we extract the corresponding transport coefficient. We reiterate, first of all, that within the NRQCD Lagrangian x H corresponds to a conserved charge (particle plus antiparticle number times the rest mass). In other words, the corresponding operator commutes with the HamiltonianĤ, [Ĥ, xĤ ] = 0. Therefore the Heisenberg operator is time independent, It turns out, however, that there are contact terms at the positions τ 1 and τ 2 when considering the correlator in eq. (3.5). For τ 1 > τ 2 , which will turn out to be the case relevant for us (see below), we need to have τ > τ 1 or τ < τ 2 , otherwise the correlator vanishes up to exponentially small corrections. Therefore τ can only be chosen from the range τ ∈ [0, τ 2 ) ∪ (τ 1 , β]. In order to re-express eq. (3.5) in the canonical formalism in this case, we recall that the Euclidean path integral corresponds to a time-ordered expectation value. Now, within the NRQCD theory, the operatorsθ † andθ are the creation and annihilation operators for the quark states, respectively. For the antiquark states,χ corresponds to a creation operator andχ † to an annihilation operator. Thereforeθ †χ creates a quark-antiquark state, and χ †θ annihilates it. According to eq. (2.5), xĤ measures the rest mass of the quarks and antiquarks present.
We denote energy eigenstates involving a heavy quark and a heavy antiquark by |qq, m , whereas states involving neither are denoted by |n . The corresponding energies are denoted by E m and ǫ n , respectively. Then xĤ |qq, m = 2M |qq, m whereas xĤ |n = 0. Note that the quark-antiquark states can be either bound states or scattering states.
Inserting now completeness relations 2 and making use of eq. (3.6), the numerator of eq. (3.5) can be expressed as We have restricted to the contribution from sectors of the Hilbert space with at most one quark-antiquark state present. The quark-antiquark states had to be placed at the outer boundaries in order for xĤ to give a non-zero contribution. We note that had we chosen τ 2 > τ 1 instead, then the creation operator (θ †χ )(τ 1 , 0) would have appeared to the right of the annihilation operator (χ †θ )(τ 2 , 0). Then the states |n would have to be replaced by states containing two heavy quarks and two antiquarks. This is the case irrespective of the range chosen for τ , because one of theĤ-operators is always at τ = 0. The corresponding contribution to Γ chem would be exponentially suppressed by exp(−2M/T ) and will be omitted.
Eq. (3.7) should clarify the physical meaning of eq. (3.5). As indicated by the matrix elements squared, we are considering the decays of quark-antiquark states into light degrees of freedom. Physically, the quark-antiquark states are assumed to be initially close to thermal equilibrium; our general framework corresponds to a linear response analysis, so that in eq. (3.7) the quark-antiquark states are in thermal equilibrium.
Understanding the time dependence of eq. (3.7) requires care. Let us denote C mn ≡ 4M 2 Z e −βEm qq, m|θ †χ |n n|χ †θ |qq, m . Recalling the restrictions on τ the function in eq. (3.5) can then be written as Taking a time derivative; carrying out one integration with the help of the resulting Dirac δ-function; and representing the denominator as 1/( This can be integrated into where the (infinite) integration constant can be omitted, because its Fourier transform has no non-trivial cut. Going over to the normalization of eq. (3.4), a Fourier transform according to eq. (2.10) yields (3.12) The corresponding spectral function reads where we noted that the dynamical energy scales contained in the effective description, reflected by the energy eigenvalues ǫ n , are by assumption much smaller than the heavy quarkantiquark energies E m ∼ 2M , so that the contribution emerges from s ∼ 2M . 3 Taking finally 3 This is a subtle point. It could be said that the energies Em ∼ 2M had already been integrated out in order to arrive at the NRQCD description. However, in order to get the correct limit, we need to keep the rest mass in the heavy quark Lagrangian here (cf. eq. (A.3)); it can only be shown a posteriori (see below) that the rest mass can be shifted away also in our finite-temperature observables.
the limit in eq. (2.11) and inserting subsequently the value of C mn from eq. (3.8) yields Here an essential point was that there was no functional dependence on the energy eigenvalues ǫ n , so that we could re-identify the sum n |n n| as a unit operator. Subsequently we returned to a path integral representation. We also indicated one of the time arguments by 0 + in order to maintain the correct ordering. Eq. (3.14) is one of our main results. Effectively, it represents the expectation value of the imaginary part of the singlet operator in eq. (2.7), albeit in a spacetime with a Euclidean signature, with a periodic time direction, and with a specific time ordering. Eq. (3.14) shows that even though the general representation of the absorptive parts of 4-quark operators is cumbersome (cf. sec. 3.1), the final transport coefficient can be extracted from an "almost" local 4-fermion operator. The time ordering in eq. (3.14), represented by the time argument 0 + , has however a specific consequence, to which we now turn.

Wick contractions
In order to express eq. (3.14) in terms of propagators, we need to recall some of their basic properties. The heavy quark propagator is defined in appendix A. An important consequence of the fact that non-relativistic propagators are defined by equations which are of first order in time is that G θ (τ 2 , x; τ 1 , y) = θ(τ 2 , x)θ † (τ 1 , y) is obtained by integrating into the region τ 2 > τ 1 (cf. eq. (A.3)), and G χ (τ 2 , x; τ 1 , y) = χ(τ 2 , x)χ † (τ 1 , y) is obtained by integrating into the region τ 1 > τ 2 . However in eq. (3.14) the time argument of θ † is larger than that of θ, and the time argument of χ is larger than that of χ † . Therefore a non-zero contraction requires integrating around the imaginary-time direction. The situation is illustrated in fig. 1.
With these observations and the relations in eqs. (A.5) and (A.6) in mind, the correlator in eq. (3.14) can be expressed as where α, γ ∈ {1, ..., N c } are colour indices, and i, j ∈ {1, 2} are spin indices. We may also express χ f , entering eq. (2.16), in the same notation. A few manipulations, including the use of eq. (A.6), lead to Inserting eq. (3.15) into eq. (3.14) and normalizing the result according to eq. (2.16), with χ f inserted from (3.16), we get our final result for the non-perturbative Sommerfeld factor: (3.17) The factor 2N c corresponds to the dimension of the propagator matrices.
The following physical interpretation can be suggested. Effectively we are computing the thermal expectation value of the imaginary part of the 4-quark operator, cf. eq. (3.14). But the contractions are non-trivial in that a quark-antiquark pair is generated at time 0 and annihilated at time β. The circling of the imaginary-time direction guarantees that the corresponding physical states are thermalized. The system is allowed to decide whether the propagation takes place in the form of open or bound states, as long as they are distributed according to the proper thermal weights. This physics is normalized by the propagators of two independent heavy quarks or antiquarks.
In order to write down the corresponding correlator for the octet case, eq. (2.9), we define the ratios Summing over the generators T a and normalizing the result according to eq. (2.19) we get (3.21) With the same notation,S 1 = P 2 /P 2 1 . Note that all the quantities entering the singlet and octet factors are manifestly gauge-invariant.

Perturbative estimates
In order to evaluate eq. (3.17) within (resummed) perturbation theory, we rewrite it in the form of a spectral representation. Denoting by ρ qq (E, k) the spectral density of quarkantiquark states, and representing the spectral density of the single-particle states in the denominator by the non-relativistic free form ρ q (E, p) = 2N c πδ(E − E p ), the spectral representation of eq. (3.17) reads 4S .

(4.1)
Here k is the momentum of the quark-antiquark pair with respect to the heat bath. The Boltzmann weight e −βE implies that states with a small mass appear with a high likelihood. We evaluate eq. (4.1) in three ways, with increasing sophistication. At tree-level, the spectral function originates from free scattering states, and reads Inserting this into eq. (4.1) and carrying out the integrals over E and k, immediately yields S 1 = 1. For future reference, it is helpful to express the tree-level result in another way as well. Making use of the non-relativistic form E p = M + p 2 /(2M ), the integrals over p, q in eq. (4.2) are readily carried out, yielding  [26,27]) pseudoscalar spectral function in hot QCD, with E ′ ≡ E − 2M − k 2 /(4M ) denoting the energy with respect to the 2-particle threshold. Here k = |k| is the spatial momentum with respect to the heat bath. The tree-level result multiplied by the Sommerfeld factor S 1 from eq. (4.6) reproduces the above-threshold spectral shape, but it misses the (thermally broadened) resonance-like contribution below the threshold.
This result is illustrated with a dashed line in fig. 2. We note that it amounts to −1/3 times the vector channel spectral function studied in ref. [26], i.e. the pseudoscalar channel [27]. As a second step, we include the Sommerfeld enhancement in its standard form. Employing QCD language, we refer to this approximation as a contribution from "open" (or "abovethreshold", or "scattering") states. Then The Sommerfeld factor reads [3] We note in passing that expressing both the numerator and denominator of eq. (4.1) in center-of-mass variables, and inserting ρ open from eq. (4.4), a simple computation yields This agrees with the expression given in eq. (2.15). As a third level of sophistication, we evaluate eq. (4.1) by including the full spectral shape shown in fig. 2. We refer to this as a contribution from "open+bound" states. Substituting variables as E ′ = E − 2M − k 2 /(4M ); noting that the 1-loop resummed spectral function 5 depends on E ′ only [26], so that we can write ρ qq (E, k) = ρ(E ′ ); and carrying out both integrals in the denominator of eq. (4.1) as well as the integral over k in the numerator, 5 The 1-loop resummed spectral function corresponds to the imaginary part of a Coulomb Green's function, with the potential appearing in the inhomogoneous Schrödinger equation computed up to order g 2 in Hard Thermal Loop resummed perturbation theory. This potential has a real part, accounting for "virtual" thermal corrections such as Debye screening, as well as an imaginary part, accounting for "real" thermal corrections such as elastic 2 → 2 scatterings with light plasma constituents, which lead to thermal broadening.
The cutoff scale Λ ′ has no practical significance, given that the spectral function vanishes for E ′ ≪ −α 2 s M (cf. fig. 2). A numerical evaluation of eq. (4.8), based on the spectral function shown in fig. 2, is illustrated in fig. 3. It is obvious that at low temperatures, the bound state contribution completely dominates the result, boostingS 1 by more than two orders of magnitude.
In the octet case, where the interaction is repulsive, we only include the contribution from the open states, like in eq. (4.7). The octet Sommerfeld factor reads [3] Numerical results are displayed in fig. 3, showing a modest suppression.

Lattice analysis
We have measured the observables defined in eqs. (3.17)-(3.21) using lattice NRQCD. The gauge configuration ensemble corresponds to N f = 2 + 1 dynamical quark flavors at a temperature of a few hundred MeV. In this exploratory study, a single lattice spacing and a single spatial volume were employed. The same thermal gauge configurations have previously been used for studying in-medium modifications of bottomonium spectral functions [31] and the electric conductivity of the quark-gluon plasma [32]. The configurations used correspond to an anisotropic lattice setup, with the spatial lattice spacing a s and the temporal one a t related through a s /a t ≈ 3.5. The tuning of the anisotropy parameter and the vacuum properties of the lattice action were studied by the Hadron Spectrum Collaboration [33,34], yielding a s = 0.1227(8) fm, M π ≃ 400 MeV, M K ≃ 500 MeV. The critical temperature was estimated from the behaviour of the Polyakov loop expectation value [35], with the result T c = 185(4) MeV. The heavy quark mass for the physical bottom case, a s M kin = 2.92, was tuned by extracting a kinetic mass from dispersion relations and matching it onto the spin-averaged mass of the η b and Υ mesons. Further details on the lattice setup can be found in ref. [31].
It is important to note that the contribution of the bare heavy quark "rest mass" drops out in the ratios in eqs. (3.17) and (3.21), and radiatively generated self-energy divergences drop out as well. Therefore we can leave out the bare rest mass from the equation of motion defining the heavy quark propagator, which facilitates the measurement. The form of the NRQCD propagator is specified in appendix A.  the lattice spacings are asymmetric, with a s /a τ ≈ 3.5. The same configurations (10 3 per parameter set) were previously used e.g. in refs. [31,32]. Errors are statistical only; forS 1 andS 8 they were obtained with a jackknife analysis. The results of this table are illustrated in fig. 4.
The main results of our lattice study are given in table 1 and in fig. 4. We note that in addition to our main N f = 2+1 lattice calculation, we have also carried out tests on quenched lattice gauge configurations with a symmetric lattice spacing on a small lattice volume 4 × 8 3 . For completeness the result of quenched test is shown in table 2.
Comparing our N f = 2 + 1 lattice results in fig. 4 with the perturbative results in fig. 3, the following conclusions can be drawn: (i) The lattice and perturbative results forS 8 are of similar magnitudes, and at T > ∼ 1.5T c both are somewhat below unity. At smaller T /T c the central values of the lattice results are larger, however the observable in eq. (3.21) involves a subtraction in the numerator and thus a large cancellation at T < ∼ 1.5T c . This is likely to lead to systematic uncertainties from lattice artifacts which may be much larger than the statistical ones.
(ii) The lattice results forS 1 can exceed 10 3 , whereas the perturbative ones according to the prescription used in the literature, which we have denoted byS 1,open (cf. eq. (4.7)), never exceed 10. Even though the perturbative result may be expected to have ∼ 50% uncertainties (see below), such a discrepancy of two orders of magnitude cannot be accounted for by conceivable perturbative uncertainties.
(iii) The lattice results forS 1 are in qualitative agreement with the perturbative ones in- T /T c a s M kin M kin /T 10 2 P 1 10 3 P 2 10 4 P 3S1S8  cluding bound state contributions, which we have denoted byS 1,open+bound . Certain differences do remain: the mass dependence at fixed T is more rapid in the perturbative results (cf. right panels of figs. 3 and 4), and the temperature scales are somewhat shifted (cf. left panels of figs. 3 and 4). We stress, however, that the vertical axes of the plots are logarithmic; without the inclusion of bound state contributions, there would be differences of up to two orders of magnitude. With the inclusion of bound states, orders of magnitude can be accounted for. For the remaining discrepancies there are many likely explanations. First of all our leading-order perturbative computation has uncertainties on the ∼ 50% level. Amongst others, the perturbative predictions are sensitive to the parameter values used, in particular the scale chosen for the gauge coupling and the scheme chosen for the quark mass; in the absence of an actual higherorder computation ofS 1 , the choices made amount to ad hoc recipes. Second, the lattice results correspond to quark masses larger than in the continuum computation, which implies that physical scales cannot be compared unambiguously (this is reflected by the different values of T c /MeV). Third, the lattice results have no continuum limit and take place in a finite volume, V ≈ (2.9 fm) 3 . It would be interesting to carry out a refined lattice analysis and a systematic higher-order perturbative computation in order to see if the discrepancies get thereby reduced.

Discussion and outlook
We have presented a power-counting argument showing that bound states (if they exist) dominate the thermal co-annihilation rate of kinetically equilibrated non-relativistic particles at low temperatures (discussion around eq. (2.3)), as well as a theoretical derivation of basic formulae (eqs. Even though our study was inspired by the co-annihilation rate of weakly interacting nonrelativistic particles in cosmology, which determines their freeze-out temperature, we applied the formalism to obtain a non-perturbative estimate of the heavy quark chemical equilibration rate in QCD, relevant for heavy ion collision experiments. Our numerical investigation confirmed the existence of Sommerfeld enhancement in the case of attractive interactions as well as a reduction (at T > ∼ 1.5T c ) in the case of repulsive interactions (cf. fig. 4). At T > ∼ 1.5T c the repulsive case shows reasonable agreement with perturbation theory (cf. fig. 3).
In the attractive case, corresponding to the factorS 1 , we found a much larger Sommerfeld enhancement than predicted by the perturbative formula sometimes used in literature (cf. eq. (4.7)). In accordance with the power-counting argument around eq. (2.3) we believe that the difference originates from bound state contributions, which were omitted in eq. (4.7). 6 Due to their smaller mass, bound states appear with a less suppressed Boltzmann weight in the thermal ensemble than scattering states. Note also that bound state formation is a relatively speaking fast process, without any exponential suppression factors; thermal friction, caused by 2 → 2 scatterings with light plasma constituents, lets an open qq pair lose energy until it appears with its proper thermal weight. Including bound state decays through a numerical determination of the imaginary part of a thermal Green's function, ρ qq (E, k), we did find qualitative agreement between perturbative and lattice results (cf. figs. 3 and 4). 2. This means that we have been deeper in the enhancement regime than is typical for cosmology. That said, there are cosmological models in which bound state formation has been demonstrated in vacuum (cf. e.g. refs. [37][38][39][40]). It might be interesting to use our formalism to study whether or not bound states could affect the thermal freeze-out process in these cases. Another direction worth a look are models based on strongly interacting dark sectors (cf. e.g. refs. [41][42][43][44] and references therein).
Let us finally return from cosmology to hot QCD. Inserting the values of the absorptive parts of the 4-quark coefficients from ref. [15] and assuming that the octet Sommerfeld factors are spin-independent, the heavy quark chemical equilibration rate evaluates to [13] Γ chem ≃ Inserting the highest temperature reached in current heavy-ion collisions, T ∼ 400 MeV, a charm quark mass M ∼ 1.5 GeV, a running coupling α s ∼ 0.25, and valuesS 1 ∼ 15, S 8 ∼ 0.8 as estimated in the current study for N f = 3, we get Γ −1 chem ∼ 150 fm/c. Therefore heavy quark chemical equilibration is unlikely to take place within the lifetime ∼ 10 fm/c of the current generation of heavy ion collision experiments. This can be compared with their kinetic equilibration time scale which could be as small as ∼ 1 fm/c in the case of charm quarks [45,46]. It may be noted, however, that Γ chem (eq. (6.1)) changes rapidly with temperature. Already a modest increase, in the ballpark of the planned Future Circular Collider (FCC) heavy ion program [47], could therefore yield chemically equilibrated charm quarks.
On the lattice the non-relativistic quark propagator is calculated from (· ≡ 0, 0) where U 0 is a time-direction gauge link. The lowest-order Hamiltonian reads where ∆ (2) is a discretized gauge Laplacian. The higher order correction is where g 0 is the bare gauge coupling, and unspecified notation is explained in ref. [31]. The parameter n, which stabilizes the high-momentum behaviour of the propagator, is set to 1 for a s M kin = 2.92 and to 3 for a s M kin = 1.5. The electric and magnetic fields appearing in eq. (A.11) were implemented in a tadpole-improved [48] form, with the improvement factor u s = 0.7336 for the spatial link and u τ = 1.0 for the time link [33].