Debye screening mass of hot Yang-Mills theory to three-loop order

Building upon our earlier work, we compute a Debye mass of finite-temperature Yang-Mills theory to three-loop order. As an application, we determine a $g^7$ contribution to the thermodynamic pressure of hot QCD.


Introduction
In electromagnetic plasmas, the Debye screening mass m el -or the inverse screening length of electric fields within the plasma -is a most fundamental quantity. While it is sometimes defined as the small-momentum limit of the static Coulomb propagator 1/[k 2 +Π 00 (0, k)] as m 2 D = Π 00 (ω = 0, k → 0), where Π 00 (ω, k) is the longitudinal part of the photon self-energy, an alternative definition is in terms of the pole of the same static propagator, (1.1) Equivalently, in a quark-gluon plasma, the Debye screening mass parameterizes the dynamically generated screening of chromo-electric fields, due to the strong interactions of hot quantum chromodynamics (QCD). Now, since the longitudinal gluon self-energy is not a gauge invariant quantity, and since its static low-momentum limit exhibits severe infrared (IR) divergences [1,2], it is not at all obvious that the above definitions for a physical quantity make sense. Indeed, after a number of investigations of this matter [3][4][5][6][7], it has JHEP11(2015)121 turned out that the definition of eq. (1.1) is the physically sensible one, 1 leading to a gauge invariant and infrared finite Debye mass also at higher orders in perturbation theory, as has been demonstrated in refs. [9][10][11][12]. The analytic treatment of such hot QCD systems can be transparently organized after identifying the different dynamically generated energy scales. Applying the concept of effective theories (EFTs) to this multi-scale system, it has been understood how to reduce the root cause of the IR problem to a well-defined non-perturbative lattice measurement which, after mapping to the continuum, constitutes a systematic approach that evades the IR problem and renders thermodynamic observables theoretically computable, allowing for systematic improvements. In the case of QCD in thermal equilibrium, one can identify three relevant energy scales, being effectively described by a set of three distinct theories: 4-dimensional hot QCD, 3-dimensional Electrostatic QCD (EQCD) and 3dimensional Magnetostatic QCD (MQCD), respectively. Together, after proper matching of their parameters, they allow for consistent weak-coupling expansions of static quantities.
In the present paper, by a three-loop determination of the Debye mass of hot Yang-Mills theory, we intend to contribute yet another coefficient to the weak-coupling EFT setup, which is needed for matching QCD and EQCD. Our motivation is threefold. First, our new result immediately determines a (gauge invariant part of the) g 7 contribution to the thermodynamic pressure. The importance of this higher-order contribution lies in the fact that it represents the next-to leading order (NLO) correction to what can be called the physical leading order, 2 ultimately enabling the first sound statements about convergence and (renormalization) scale dependence, going beyond existing discussions that are based on truncated (or, in the modern understanding, incomplete 3 LO) versions of the series.
Second, the Debye mass plays a prominent role in various channels of gauge-invariant gluonic screening masses [7,14,15]. Within a weak-coupling expansion, a number of these are at leading order proportional to the Debye mass followed by a formally sub-leading, but numerically large, logarithmic correction. There are examples, however, where the logarithmic term is known to be small, such that the Debye mass dominates the functional behavior in the phenomenologically relevant temperature range, such as has been observed for the lowest-lying color-magnetic screening mass [15].
Third, we regard the determination of the 3-loop Debye screening mass m 2 E , as an EFT matching parameter, as a proof-of-principle that also the corresponding evaluation of the 3-loop effective gauge coupling constant g 2 E is within reach. The latter does bring a direct phenomenological application, allowing for a precise comparison of the so-called spatial string tension σ s (as evaluated in the EFT setting) with lattice determinations, where a previous 2-loop analysis had shown great promise, while leaving room for further corrections. Furthermore, such a comparison can be regarded as an important consistency check, validating the EFT approach as a whole. 1 Note that the proper definition of the pole location is a bit more subtle, see e.g. section 5 of [8], but this makes no difference here. 2 In the EFT setup, this is defined to be the order in which all potentially large logarithmic terms originating from the different physical scales of the problem have entered the observable at hand for the first time; for the pressure, this happens to be the order g 6 [13], see also section 4 below. 3 Since the effect of large logarithms is then not considered properly.

JHEP11(2015)121
The structure of the paper is the following. In section 2 we provide a brief overview of the theoretical setup and of the matching relations between full QCD and EQCD and sketch the status of the Debye mass computation reached in ref. [16]. In section 3 we complete the calculation, restricting ourselves to the case N f = 0 (pure gauge theory), express the bare mass parameter in terms of a few master sum-integrals, renormalize and analyze our result numerically. We then use the renormalized result for extracting a g 7 contribution to the QCD pressure in section 4. We conclude in section 5, while some technicalities are relegated to the appendices.

Setup: effective theory and matching
At finite temperatures, gluons exhibit three characteristic momentum scales (2πT , gT and g 2 T ) of which the ultra-soft color-magnetic mode (g 2 T ) leads to a breakdown of the ordinary perturbative expansion [2]. A straightforward approach that sidesteps this problem is scale separation. This is achieved by constructing dimensionally reduced effective theories [17][18][19] whose parameters are matched to full QCD. Following this procedure, two effective theories in d = 3 dimensions emerge: electrostatic QCD (EQCD) is an SU(N ) gauge theory coupled to a massive scalar field in the adjoint representation. EQCD contains two scales (the soft and the ultra-soft one) while the hard scale enters only through the perturbative matching of the theory's parameters. Magnetostatic QCD (MQCD) is a pure SU(N ) gauge theory containing only the non-perturbative ultra-soft scale, whereas the hard and the soft scales enter again only through the matched parameters of the theory. MQCD can only be studied with non-perturbative methods such as lattice QCD [20,21].
In the following, we are concerned only with matching hot QCD to EQCD. The bare EQCD Lagrangian reads The operators quartic in the fields A 0 are linearly independent for N > 3 only, while for N = 2, 3 we have TrA 4 0 = 1 2 (TrA 2 0 ) 2 . The term δL E collects the infinite tower of higherorder operators that are generated by integrating out the hard scale, the lowest of which have been classified in ref. [22]. We shall be needing a single one of them, cf. eq. (2.8) below.
The detailed framework of performing the matching computation has been presented in [16,23]. Here, we merely provide a concise version of it and generalize the matching condition in order to account for higher-order operators. The general prescription is to require that various static quantities computed in both theories match to a certain order in a strict perturbative expansion with respect to the gauge coupling g. By using the background field gauge, we make sure that on the QCD side only the coupling constant g requires renormalization [24,25].

Screening mass in QCD
Following eq. (1.1), we define the screening mass m el as the pole of the static (K 0 = 0) momentum-space propagator of A 0 with the on-shell condition k 2 = −m 2 el . On the full QCD side, writing Π ab 00 (0, k) ≡ δ ab Π E (0, k) we therefore have (2. 2) The self-energy is written as an expansion in both the gauge coupling g and in the external momentum k in order to permit a strict perturbative expansion. The k-expansion is justified due to the soft scale |k| ∝ gT at which the pole in the propagator arises [26]: Solving iteratively, eq. (2.2) leads to the following expression for the screening mass: Full d-dimensional representations of the various coefficients Π En can be found in appendix C of [16]. For the reader's convenience, we have collected the corresponding oneand two-loop self-energies in appendix 6 below. The evaluation of the QCD self-energy tensor Π µν (K) to three-loop order generates approximately 500 Feynman diagrams, necessitating an automatized procedure to handle this task (cf. [16,23] and references therein). Generation of the Feynman diagrams, the color algebra computation of SU(N ), the Lorentz contraction and the Taylor expansion into external momentum have been performed using specialized software (we have employed QGRAF [27,28] and FORM [29,30]). The resulting ≈ 10 7 sum-integrals have then been reduced via systematic use of integration-by-parts (IBP) relations [31] in the thermal context [32] to a small number of master sum-integrals J and I, as defined in appendix 7.
For N f = 0, the bare 3-loop result of refs. [16,33] reads with pre-factors α 1...10 given in appendix 8. Inspecting this result, there are two trivial products of one-loop tadpole sum-integrals ∼ I · I · I that are already known analytically (cf. eq. (7.3)) and eight non-trivial three-loop cases J αβγ ab00cd of basketball-type, of which so far only a few terms of the one multiplying α 1 (originally given in [34] and re-evaluated as a specific case of the class J 000 N 10011 in [35]) as well as the one multiplying α 4 (treated as a special case of J M 00 N 10011 in section 5 of [16]) are known. As discussed in ref. [16], the difficulty in calculating the 3-loop sum-integrals J i that appear in eq. (2.5) lies in the fact they would have to be expanded beyond the constant JHEP11(2015)121 term in ǫ (in our notation, d = 3 − 2ǫ) due to their singular pre-factors, cf. appendix 8. Conventional techniques of evaluating basketball-type sum-integrals [36,37], relying on setting d = 3 for determining their constant parts in coordinate space representations, make this task difficult (if not impossible). We will show in section 3 how to proceed.

Screening mass in EQCD
In EQCD, writing the self-energy of the adjoint scalar A 0 as Π ab EQCD = δ ab Π EQCD , the pole mass of the A 0 propagator is By performing the same twofold expansion of the A 0 self-energy Π EQCD as in eq. (2.3), it vanishes to all orders in perturbation theory due to the absence of any scale in the resulting three-dimensional vacuum integrals 4 (cf. discussion in section 2.1 of [16]). However, in [16] higher-order operators to eq. (2.1), such as those computed in [22] and that contribute to were not yet considered. Indeed, in the context of the double expansion, in which only scale-free vacuum integrals emerge, one might expect that these higher-order contributions all vanish.
It turns out, however, that one higher-order operator does contribute to the self-energy, since it generates a tree-level two-point contribution that is not affected by the momentum expansion. This dimension six operator can be extracted from ref. [22]: where g 2 R is the (dimensionless) renormalized 4d QCD gauge coupling. Since at the pole, the four derivatives will scale as (k 2 ) 2 ∼ g 4 T 4 , eq. (2.7) receives a correction of O(g 6 ).

Completing the calculation
As already emphasized, a technically challenging part of the matching is the evaluation of the master sum-integrals. We also mentioned the impossibility of computing sum-integrals beyond the constant term in ǫ with state of the art techniques. However, in eq. (2.5) the singular pre-factors in ǫ multiplying the master sum-integrals require their evaluation including O(ǫ).

Basis transformation
A possible way out is to search for a suitable basis transformation that -much in the spirit of the ǫ-finite basis advocated in ref. [38] -removes the singularities of the pre-factors in eq. (2.5) for the price of introducing master sum-integrals that are of a different topology and might therefore be more difficult to compute. Using the master integrals defined in appendix 7 as well as the IBP relations of appendix 9, we have succeeded in rewriting the bosonic part of Π E3 shown in eq. (2.5) as the remarkably compact expression Using the lower-order self-energies listed in appendix 6, eq. (2.4) then immediately gives Note that all dependence on the gauge parameter ξ has duly canceled in this d-dimensional result. By plugging the leading term of eq. (3.3) into eq. (2.9), we finally obtain The remaining task is to evaluate the sum-integrals J 11,12,13 that enter our expression for m 2 el and that are depicted in figure 1. To this end, the 1-loop substructure can be JHEP11(2015)121 exploited, a method pioneered by Arnold and Zhai [36,37]. Their technique of solving basketball-type and spectacle-type sum-integrals relies on a careful subtraction of subdivergences which is specific for every sum-integral in part. Nevertheless, it was possible to develop a semi-automatized procedure for an analytic calculation of the divergent parts of a large class of spectacle-type sum-integrals. This was necessary for evaluating the Mercedes type sum-integral with two inverse propagators, J 13 . Its computation required the use of the dimensional method of Tarasov [39], in which the tensor structure of a sum-integral is translated into a sum of higher dimensional scalar integrals. The last remaining pieces of the matching computation, the 3-loop vacuum sumintegrals J 12 and J 13 , have been determined only recently in refs. [40,41] and are listed in appendix 7.

Renormalization
We now turn to renormalized quantities. In the 4-dimensional theory, we need to renormalize the gauge coupling. The bare coupling g B (denoted as g in the previous sections) is related to the renormalized coupling g R via g 2 B = g 2 R µ 2ǫ Z g (note that g 2 R is dimensionless), µ being the MS scheme scale defined asμ 2 = 4πe −γ E µ 2 and where, for N f = 0, In the 3-dimensional theory, only the mass parameter requires renormalization. This is due to the fact that the EQCD Lagrangian in eq. (2.1) is super-renormalizable in d = 3−2ǫ dimensions. The few divergent diagrams in this theory arise solely at two-loop order and account only for the mass renormalization. Hence, the renormalization relations are simply the renormalized EQCD couplings g 2 ER and λ (1/2) ER being of dimension one, while the EQCD mass counterterm reads [42][43][44] This leads to an exact RG equation for the mass parameter: µ 2 3 ∂ µ 2 3 m 2 ER = 2ǫ δm 2 E . In order to express δm 2 E in terms of the 4d QCD coupling g 2 R , the EQCD parameters g 2 E and λ (1/2) E need to be matched to full QCD at leading order (1-loop). The corresponding relations are given in the literature as [18,45], however without regularization parameter ǫ, which we here need due to the presence of the 1/ǫ divergence in eq. (3.7). We hence need d-dimensional JHEP11(2015)121 matching relations, which can be extracted from [19], in our notation, as Using these expressions, the mass counterterm simplifies to such that we finally get the renormalized Debye mass up to 3-loop order as where the 1-loop and 2-loop contributions have been calculated in ref. [46] and where we have expressed the contribution coming from eq. (2.8) separately. The g 6 R coefficient contains two independent mass scales, µ and µ 3 . The first arises through the 4d dimensional regularization scheme and the subsequent renormalization of the 4d coupling, g R , whereas the second scale µ 3 enters through the regularization of the divergent integrals of the 3d SU(N ) + adjoint Higgs theory and ultimately through the mass renormalization.

Numerical evaluation
For the numerical evaluation of m ER , the running of the 4d coupling with respect to the energy scale is obtained by solving the RGE iteratively to three-loop order [47,48] with t = ln[μ 2 /Λ 2 MS ] and Λ MS the QCD scale defined in the MS scheme [48,49]. In order to display numerical results, we need to choose values for the two arbitrary mass scales, µ and µ 3 . For the former one, we adopt the procedure of of minimal sensitivity [46]. The scale is computed to beμ opt /T = 4πe −γ E − 1 22 ≈ 2π. We then extend this idea to choosing µ 3,opt . As m 2 ER (μ, µ 3 ) is not a monotonic function with respect toμ, we impose that the absolute variation of m 2 ER (μ, µ 3 ) in the interval (μ opt /2, 2μ opt ) is minimal for a specific scale µ 3 ≡ µ 3,opt , or Taking the absolute variation ensures that an oscillatory behavior of m 2 ER in the considered interval is ruled out to be regarded as a scale for minimal sensitivity. Solving the equation numerically, we obtain µ 3,opt /T ≈ 2.85. In figure 2, we analyze the running of the Debye mass with respect to the temperature (T ). We have used the solution of the renormalization group equation of the 4d coupling g, in which the QCD β-function was truncated after O(g 8 ). The parameter Λ MS corresponds to the QCD scale defined in ref. [46], and here simply sets the scale through the t-dependence of the running coupling eq. (3.11). The arbitrary scaleμ was chosen at the point where the effective coupling g ER has a minimal sensitivity to it:μ opt /T ≈ 2π. The 3 bands in the plot arise by varyingμ = (0.5 . . . 2) ×μ opt . Using the prescription described above for choosing a sensible scale for the EQCD scale parameter, µ 3 /T ≈ 2.85, we obtain a three-loop result with a vanishingly narrow band width.

JHEP11(2015)121
From the figure, one notices a slight increase of the Debye mass with respect to the 2-loop result. In addition, the sensitivity with respect to the arbitrary scaleμ decreases, which indicates that the perturbative expansion up to 3-loop order shows good convergence properties.

A g 7 contribution to the QCD pressure
Having the 3-loop Debye mass at hand, we can use it to extract a gauge-invariant piece of a higher-order perturbative (g 7 ) correction to the QCD pressure. The order g 7 owes its importance to the fact that it represents, in the effective theory setup we are working in, the leading correction to what has been called the physical leading order, i.e. all terms up to order g 6 that, for the first time, include all potentially large logarithms entering the QCD pressure [13].
As already discussed in the introduction, in the effective theory framework, the QCD partition function factorizes, such that the pressure p QCD splits into three parts p E , p M and p G originating from contributions of the hard-2πT , soft-gT and ultra-soft scale g 2 T , respectively [26,50]. These three parts can be extracted as matching coefficients of QCD,

JHEP11(2015)121
EQCD and MQCD, along the following chain of equations (recall d = 3 − 2ǫ), supplemented with the corresponding matching of the couplings, as has been explained in the previous sections on the example of the chromo-electric screening mass. It turns out that the different pieces in eq. (4.1), when re-expressed in terms of the renormalized 4d gauge coupling g R (and omitting logarithms of the coupling), contribute as Both parts of the pressure coming from the 3d reduced effective theories describing soft and ultra-soft effects, EQCD and MQCD, contain contributions of order g 7 R . The latter part originates from matching of the MQCD gauge coupling g M , as p G = T g 6 M c G , with c G being a non-perturbative coefficient determined in [13,51], and the coupling g 2 M ∼ g 2 ) with coefficients known from e.g. [44]. The contributions of order g 7 R to the pressure of hot QCD from the soft scale gT entering through EQCD originate from a number of different sources. Within EQCD, the expansion reads (see e.g. [26]) with known 2. . . 4-loop coefficients a 2...4 [52], and a 5 coming from a 5-loop calculation of the EQCD pressure including contributions from higher-order operators [22] omitted in eq. (2.1), both of which remain unknown to date. The expansion parameters above relate to the 4d coupling as such that we can already fix the other g 7 pieces by multiplying out the expansion parameters and our new 3-loop result for m 2 E . Writing

JHEP11(2015)121
where α E8 represents our new 3-loop result of eq. (3.10), this coefficient contributes a g 7 term to the QCD pressure, via with logarithms L and L 3 as in eq. (3.10) above.

Conclusions
In this paper, we have determined the Debye screening mass of hot Yang-Mills theory to 3-loop (NNLO) accuracy, by combining various results of a long-term project, with a number of independent ingredients, each of which needed novel state-of-the art techniques for a successful determination. While the feasibility of such a precision calculation in the thermal context was not at all clear from the outset, we have here succeeded to overcome the last major obstacle -namely to map a sum of seven non-trivial master sum-integrals (that remained after the IBP reduction algorithm had halted, but which would have been needed to an expansion depth for which no technique existed) to a small number of (three) computable cases. We have achieved this reduction in complexity of the calculation by a clever basis transformation, which was made possible by searching our extensive database of IBP relations. To our utmost satisfaction, the final assembly of all building blocks revealed (a) gauge parameter independence, (b) a finite result after renormalization, and (c) good convergence properties. We were therefore able to add another term to the pool of known (and heavily used) matching coefficients of EQCD, and to utilize this fresh term right away, determining one of the physical next-to-leading order (g 7 ) contributions to the pressure of hot QCD. As we have discussed above, this is of course not the complete g 7 result, but represents a well-defined (and gauge invariant) contribution to it.
Thus, looking back on the many technical and systematic advances that have been made during this project, we conclude that a determination of the 3-loop effective gauge coupling g 2 E , which originates from a (by one) higher moment of two-point functions and which is hence amenable to the same techniques as m 2 E , should be within reach. Another avenue for further investigations would be a generalization of our strategy to fermionic contributions, which we have ignored completely here, setting N f = 0. While the reduction to basketball-type master integrals is done [16], open problems are finding a suitable basis change, and evaluating the corresponding fermionic masters -which, containing no zero-modes on fermionic lines, could however turn out to be less involved than the bosonic cases that we have used here.

JHEP11(2015)121
The result for the Debye mass shows a good convergence in a large temperature range, suggesting that already at the three-loop order the corrections are numerically negligible, but would serve merely in future perturbative calculations (of observables such as the pressure of QCD) to ensure finite renormalized results (cancellation of UV divergences). In the light of the apparent fast convergence of the analytic result, a re-evaluation of the non-perturbative constant of the QCD screening mass can be considered [53,54] since this latter quantity is in fact used in studies of quark gluon plasma parameters such as the jet quenching parameter [55,56].

Master sum-integrals
Let us define a generic notation for massless 3-loop vacuum sum-integrals where all momenta are understood bosonic. Let us remark that from the outset, only integrals J 000 abcdef enter the calculation; however, due to the fact that the IBP relations act

JHEP11(2015)121
in the d spatial dimensions only and hence explicitly break (d+1)-dimensional rotational invariance introducing the 4-vector U = (1, 0), the numerator structure of eq. (7.1) occurs naturally in the reduction step. The original integral reduction eq. (2.5) contains only basketball-type sum-integrals J αβγ ab00ef as well as trivial products of one-loop cases Of these, we need the products containing the numbers γ E and γ 1 arising from the expansion of the Riemann Zeta function around its pole at unity, ζ(1 − ǫ) ≈ −1/ǫ + γ E + γ 1 ǫ + . . . , as well as Z n ≡ ζ ′ (−n)/ζ(−n). The specific cases on the left-hand sides of eqs. (9.1)-(9.3) are more complicated integrals, which have however already been evaluated up to their constant terms in a number of tour de force computations documented in refs. [40,58,59] and [41], respectively, from where we collect the results for convenience: 5 +γ E +2Z 1 ǫ +c 3 ǫ 2 +O(ǫ 3 ) . (7.8) The constant parts are known numerically only, and have been determined in the above-mentioned references to be c 1 ≈ +43.8676(1), c 2 ≈ +93.0894417(2) and c 3 ≈ +44.629857(1). In the present computation, however, these constant parts do not contribute, since the integrals are multiplied by pre-factors ∼ ǫ, cf. eq. (3.1).

IBP relations for basis transformation
The idea of performing a basis transformation on Π E3 translates into going some steps back into its IBP reduction [16]. The goal is to search for relations that change the coefficients of the master sum-integrals (cf. the last paragraph of appendix C in [16]) in such a way

JHEP11(2015)121
as to eliminate all factors of (d − 3) in the denominator. While it is not at all clear from the outset that this can always be achieved, it happens indeed if we use the following three automatically generated IBP relations, expressed in terms of the basis of basketball-type 3-loop master integrals defined in eq. (2.6): 2)