A QCD Debye mass in a broad temperature range

The Debye mass sets a scale for the screening of static charges and the scattering of fast charges within a gauge plasma. Inspired by its potential cosmological applications, we determine a QCD Debye mass at 2-loop order in a broad temperature range (1 GeV ... 10 TeV), demonstrating how quark mass thresholds get smoothly crossed. Along the way, integration-by-parts identities pertinent to massive loops at finite temperature are illuminated.


Introduction
If two test charges are put a distance r apart within a plasma, they influence each other with a force which is weaker than the Coulomb force in vacuum, as a result of the screening caused by the light plasma particles. The potential then takes a Yukawa form, −αe −m E r /r, where the parameter m E may be called an electric or a Debye mass. In a relativistic plasma, it is of order m E ∼ gT , where T is the temperature and g is a gauge coupling.
In the present paper we focus on strong interactions, such that g is the coupling of the SU(3) gauge force. Standard applications of the QCD Debye mass can be found in the physics of heavy ion collision experiments. However, the temperatures reached there (T ≪ 1 GeV) are so low that it is questionable whether perturbative tools are viable. Here we rather take T > ∼ 1 GeV, and consider the possible role of the QCD Debye mass in cosmology.
Given that strong interactions are in thermal equilibrium in a broad temperature range, QCD does not normally play a prominent role in cosmology. However, exceptions can be envisaged. 1 For instance, it has become popular to consider scenarios in which dark matter is but the lightest among the particles of a larger dark sector. Then it is conceivable that the dark sector may also contain particles charged under QCD (cf., e.g,. ref. [1] for a review of one such framework). At high temperatures, the pair annihilation of the QCD-charged particles would be modified by Debye screening [2]. Charged particles also experience a thermal mass shift, known as the Salpeter correction in plasma physics, ∆M T ∼ −αm E /2, which can have a relatively speaking large effect if narrow degeneracies are present in the dark sector.
Another possible application concerns the decay of heavy particles, for instance righthanded neutrinos in leptogenesis. In this case it is important to know how fast the decay products (some of which could be hadronic, produced through the Higgs channel) equilibrate kinetically [3]. This requires large-angle scattering, again sensitive to Debye screening [4]. Another relevant rate, namely that of decoherence of the decay products, originates from a difference of small-angle scatterings mediated by colour-magnetic and colour-electric fields, whereby it is non-vanishing at O(αT ) thanks to m E = 0 (cf., e.g., ref. [5]).
A third application of Debye masses is that they play a role in dimensionally reduced descriptions of the electroweak phase transition [6]. In particular, the QCD Debye mass could make a noticeable appearance if some coloured scalar field is light enough to participate in the transition dynamics (cf., e.g., refs. [7,8]). This paper is organized as follows. The definition of a Debye mass is subtle beyond leading order, so we start by specifying the concept adopted in sec. 2. The main steps and methods of the computation are described in sec. 3, and results are presented in sec. 4. We conclude in sec. 5, relegating the evaluation of massive 1-loop and 2-loop sum-integrals into appendix A.

Formulation of the problem
As mentioned at the beginning of sec. 1, the leading-order definition of a Debye mass can be related to the Yukawa screening of a static potential or, equivalently, to the thermal mass that colour-electric fields obtain. In SU(N c ) gauge theory with N f massless fermions, the classic result reads [9] When we go beyond leading order, the definition of a Debye mass is no longer unique. One possibility is to define it as the inverse of a spatial correlation length related to some gauge-invariant operator [10]. This way, the Debye mass becomes non-perturbative at nextto-leading order (NLO) [11]. However, correlation lengths depend strongly on the quantum numbers of the operator chosen. There are also other non-perturbative possibilities, related e.g. to modelling the behaviour of the static potential at intermediate distances [12].
A different strategy is to define the Debye mass as a "matching coefficient" of a low-energy description, specifically of a dimensionally reduced effective theory [13,14]. There are a number of advantages with this strategy. One is that the definition is then "universal", with the same value appearing as an ingredient in the computation of many different correlation lengths [15], or even of dynamical rates [5]. Another is that as a matching coefficient m 2 E is only sensitive to the hard scales that have been integrated out, and therefore perturbative by construction. In fact, the result is known analytically up to 3-loop order in pure Yang-Mills theory [16,17], and shows remarkable convergence down to low temperatures. Fermionic effects are for this definition currently known up to 2-loop order in the massless limit [18], and up to 1-loop level in the massive case [9]. The purpose of the current study is to extend the 2-loop result for m 2 E to include massive fermions. 2 To be explicit, the action of the dimensionally reduced effective theory, often called "Electrostatic QCD" (EQCD), reads where we are employing Euclidean conventions; and A a 0 is an adjoint scalar field. In order to determine m 2 E , it is convenient to use the background field gauge [24] as a probe. 2 We note in passing that another analogous matching coefficient is the gauge coupling of the dimensionally reduced theory, denoted by g 2 E , but its determination is technically more challenging. For pure Yang-Mills theory results are available up to 3-loop level but only in somewhat incomplete numerical form [19][20][21]. Massless fermions have been included up to 2-loop level [22], mass effects up to 1-loop level [23]. We compute the temporal 2-point function with a purely spatial external momentum, Here g 2 B is the bare coupling, which is subsequently expressed in terms of the renormalized coupling g 2 . The computation within full QCD (or the Standard Model) is matched onto a computation within the effective theory, the latter also re-expanded as a perturbative series in g 2 . However, employing dimensional regularization and Taylor-expanding in external momentum, the latter computation gives a vanishing result, given that no scales appear in the propagators. Therefore, the matching coefficient is directly given by a Taylor-expanded full theory computation, after accounting for different field normalizations (or wave function corrections) within the full and effective theories, .

Main steps of the computation
The Feynman diagrams required for determining the 2-loop fermionic contributions to m 2 E are shown in fig. 1 (we do not show gluonic diagrams as our results for them agree with ref. [18]). Apart from vertices involving the strong gauge coupling g 2 , we have for illustration also included effects from the top Yukawa coupling h 2 t , even if in practice these are small. 3 The computation is carried out by employing the gauge propagator We keep ξ as a general gauge parameter (in the background field gauge, it also appears in the cubic and quartic gauge vertices [24]), verifying that it cancels exactly at the end.
After carrying out the Wick contractions, the result can be expressed in terms of the "master" sum-integrals where P ≡ (p n , p) and p n is a Matsubara frequency. The sum-integral Σ {P } signifies that p n is fermionic, i.e. p n ≡ πT (2n + 1) with n ∈ . After renaming variables, the set of 2-loop masters can then be chosen to consist of Z 20 jkl;i , Z 11 jkl;i and Z jkl;i ≡ Z 00 jkl;i . In rare cases, we add a mass m h on the bosonic line and indicate this with the index h. The mass index i is omitted in intermediate results if this can be done without a danger of confusion.
As far as 1-loop results are concerned, we need (cf. eq. (2.4)) Eq. (3.5) is relevant because Π E1 is originally expressed as a function of bare quark masses, which are expanded as m 2 In practice, Yukawa couplings h i other than h t are omitted. Similarly, the bare gauge coupling is renormalized as g 2 . The 2-loop diagrams yield products of the 1-loop masters of eq. (3.2) as well as genuine 2-loop masters defined according to eq. (3.3). All numerators can be eliminated from 1-loop masters by making use of Z r+2 The set of masters can now be reduced by making use of integration-by-parts (IBP) identities [25], generalized to finite temperature [26]. First, inspecting and taking linear combinations, leads to relations which permit us to eliminate all quadratic powers of p n , . (3.10) Second, if we choose indices leading to two independent representations of some Z 20 jkl , we can establish relations between Z 11 jkl . Considering Z 20 212 this way, we obtain from eqs. (3.9) and (3.10) the identity By using eq. (3.9) in order to eliminate Z 20 311 , eq. (3.10) to eliminate Z 20 113 and Z 20 212 , and inserting subsequently eq. (3.11), we can remove all numerators from the sum-integrals of eq. (3.7), leading to Remarkably, IBP relations also exist between masters without any numerators. In this way, we can eliminate Z 221 , Z 211 and Z 311 in favour of Z 111 , Z 112 and Z 212 . The latter set is convenient, as it turns out that Z 111 appears with zero coefficient in d dimensions, and Z 212 can be obtained from Z 112 through a mass derivative. Thereby only one irreducible master, Z 112 , remains to be determined in detail (cf. appendix A.2). 4 The relations needed, originally found via our FORM [27] implementation of Laporta-type reduction [28], read 14) Defining ∆ P ≡ P 2 + m 2 and δ P ≡ P 2 , eqs. (3.13) and (3.15) can be verified by setting s = 1 and s = 2, respectively, in the relation whereas eq. (3.14) can be established by taking a mass derivative of eq. (3.13).
Inserting eqs. (3.13)-(3.15), the C F part gets factorized, and eq. (3.12) reduces to The various functions employed in eq. (3.18) read These are gauge-independent, i.e. no ξ appears, and also finite after the insertion of the masters from appendix A, i.e. no 1/ǫ 2 or 1/ǫ appears. For Φ (7) we have approximated the result by considering the limit m h ≪ m t , as this leads to the same basis as for the pure QCD contributions. This overestimates the magnitude of Φ (7) , but given that its effect is small even then (cf. fig. 2), the approximation can be considered conservative.

Results
Our final results for the coefficients in eq. to O(ǫ 0 ). Denoting p ≡ d 3 p (2π) 3 , ω p i ≡ p 2 + m 2 i and n F (ω) ≡ 1/[exp(ω/T ) + 1], we find where m i refer to MS masses evaluated at the renormalization scaleμ. In the limit m i ≪ T the coefficients go over into and Φ (7) ≡ (4π) 2 Φ (7) /T 2 that parametrize eq. (3.18), evaluated withμ = 2πT . Right: the QCD Debye mass as a function of the temperature. They grey band originates from varying the renormalization scale in the rangeμ = (0.5...2.0) × 2πT and gives an indication of the magnitude of higher-order corrections. The "hard scale" with which m E can be compared is ∼ 2πT . The plateau-like feature centered around T ∼ 70 GeV originates from crossing the top mass threshold.
reproducing in the first four cases the expressions obtained in ref. [18]. For m i ≫ πT , all terms containing n F are exponentially suppressed. For a numerical evaluation, we set α s (m Z ) ≃ 0.118 [29], and evolve g 2 (μ) in both directions with 5-loop running [30][31][32][33], changing N f when a threshold is crossed atμ ≃ m i , and including effects from the top Yukawa up to 3-loop order [34,35]. Quark masses are likewise evolved at 5-loop level [36,37], including effects from the top Yukawa as indicated below eq. (3.6). 5 In order to account for temperature-dependent tadpole corrections ∝ T 2 , the Higgs expectation value and m i are further scaled as v T ≃ v 0 Re 1 − (T /160 GeV) 2 where v 0 ≃ 246 GeV and the crossover temperature has been adopted from ref. [38] rather than from a perturbative computation. As the top quark Yukawa coupling plays a minor role, we have resorted to 2loop running for h 2 t , with the initial condition h 2 t (m Z ) ≃ 0.95 and running taking place only forμ > m Z . The initial value of the running top mass (before applying thermal rescaling) is

Conclusions and outlook
The goal of this technical contribution, motivated by the potential cosmological applications mentioned in sec. 1, has been to estimate a QCD Debye mass, defined as a matching coefficient of the dimensionally reduced effective theory, at temperatures between 1 GeV and 10 TeV. For this purpose we have carried out a 2-loop computation, reducing the result to a small number of exponentially convergent 1-and 2-dimensional integrals, which are readily evaluated numerically.
The most non-trivial parts of our work established the IBP relations in eqs. (3.13)-(3.15) and resolved the 2-loop master sum-integral Z 112 in appendix A.2. With these ingredients, we obtain integral representations for the various functions parametrizing our result, cf. eq. (3.18), which are shown in eqs. (4.2)-(4.7) and evaluated numerically in fig. 2(left). Putting everything together and inserting the values of running Standard Model couplings, we find that quark mass thresholds are crossed smoothly enough not to be discernible by bare eye, apart from that related to the top quark, cf. fig. 2(right).
The steps of the computation have been described on a detailed level, in order to permit for the inclusion of further massive particles if present, such as of scalar fields. Hopefully these results or techniques can find use e.g. in dark matter computations involving strongly interacting co-annihilation partners, or in precision studies of the electroweak phase transition in extensions of the Standard Model.
In the fermionic case, when the mass is non-zero, no analytic expressions are available. Even if convergent sum representations in terms of modified Bessel functions can be found, in practice it is simpler to handle integral representations, such as where p ≡ d 3 p (2π) 3 and ω p i ≡ p 2 + m 2 i . Taking a mass derivative and carrying out a partial integration gives One more mass derivative yields (this time no partial integration is possible; β ≡ 1/T )

A.2. 2-loop master Z 112
Even if in the massless limit IBP identities allow to reduce Z 112 as no such factorization has been found for m = 0. The result for a fully massive Z 111 is given in ref. [41], and one might think that Z 112 could be obtained as a mass derivative thereof, however this does not work trivially as setting the third mass to zero after the derivative leads to IR divergences (linear and logarithmic). A careful consideration is thus needed for working out the reduction of Z 112 into a convergent 2-dimensional integral representation. As a first step, let us carry out the Matsubara sums. The quadratic propagator carries a fictitious mass parameter, denoted by M 2 , as an intermediate regulator. The sum-integral splits into a vacuum part, one-cut parts, and two-cut parts, with "cut" meaning that some line is put on-shell and weighted by a thermal distribution: (A.12) Here The parameter M is set to zero for the final spatial integrals, which are treated with strict dimensional regularization (this is necessary, given that the IBP identities used for reducing the result to this basis made use of the same recipe). Two of the structures in eq. (A.7) are simple to handle, namely the vacuum part and the one-cut part with a single Bose distribution: (A.14) The remaining parts are more subtle, as they are IR divergent, in a way which is not trivially handled by dimensional regularization.
Let us start by considering Z (FF) 112 , which contains a linear IR divergence but no logarithmic one. As this integral is UV finite, we may set d = 3 and carry out the angular integral, which yields Clearly this is ill-defined around p = q. In order to find a useful representation, we make use of symmetries of the integrand, re-organizing the Fermi distributions as Considering the vacuum-like integral in eq. (A.18) as an analytic function of −p 2 and taking the real part, yields Re q 1 (q 2 −p 2 ) 2 = Re 1 8π(−p 2 ) 1/2 = 0. Alternatively, if we keep the regulator M finite, δZ (FF) 112 contains a linear divergence 6 ∝ 1/M but no logarithmic or finite part of O(M 0 ). To summarize, in strict dimensional regularization we can set δZ (FF) 112 → 0. The remaining parts, Z (F) 112 and Z (FB) 112 , contain both linear and logarithmic divergences. The logarithmic divergences cancel in the sum. We find it practical to determine the sum by keeping M finite and taking M → 0 at the end, omitting again linear divergences ∝ 1/M , which are absent in strict dimensional regularization. A rather tedious analysis then yields Given the non-triviality of the steps, it is good to check that eq. (A.6) is correctly reproduced for m/T → 0. The individual parts contain coefficients ∝ 1/m 2 , so we need to expand to O(m 2 ). The integral appearing in eq. (A. 19) can be expanded as p n F (ω p ) ω p ln me γ E 4πT + 1 + ω p 2p ln ω p + p ω p − p = T 2 24 2 + (ln ζ 2 ) ′ − ln π + 2m 2 (4π) 2 ln 2 me γ E 4πT + (1 + 2 ln 2) ln me γ E 4πT + 3 ln 2 − whereas the contribution from eq. (A.17) can be numerically verified to behave as [n F (ω p ) + n F (ω q )] 2 (ω p + ω q ) 2 − [n F (ω p ) − n F (ω q )] 2 (ω p − ω q ) 2 .