Reference : A. Bazavov, N. Brambilla, P. Petreczky, A. Vairo, J. H. Weber (TUMQCD Collaboration), /Color screening in (2+1)-flavor QCD/, Phys. Rev. D 98, 054511 (2018).

They measure, on the lattice, the correlation function of two static (infinitely heavy) colour charges separated by a distance \(r\) in a thermal bath of gluons and \(2+1\) dynamical quark flavours, over \(140 \le T \le 5814\) MeV. From that correlator they extract a free energy and ask: at what distance does the medium start to screen the colour force, and can weak-coupling effective field theories describe it?

The answer: screening switches on at \(rT \approx 0.3\), and for \(0.3 \lesssim rT \lesssim 0.6\) weak-coupling EFTs (pNRQCD, EQCD) describe the data.

word	meaning
static quark	a quark with \(m \to \infty\): it does not move, it only sits at a point and carries colour. Its worldline is a straight line in the time direction.
Euclidean time \(\tau\)	at finite temperature the time direction is imaginary and periodic with period \(\beta = 1/T\). Every field must come back to itself after \(\tau \to \tau + 1/T\).
Wilson line	the path-ordered exponential \(\mathcal{P}\exp\left(ig\int A_\mu dx^\mu\right)\). It is the object that parallel-transports colour from one point to another so that a quark-antiquark pair separated in space is a gauge-covariant object.
Polyakov loop	a Wilson line that runs all the way around the periodic time direction. Because of periodicity it closes on itself, so its trace is gauge invariant. Physically it is the worldline of one static quark living for a whole Euclidean period, i.e. "insert one static quark into the bath".
\(Z(N_c)\) centre symmetry	in pure gauge theory the action is invariant under multiplying all time-links in one time-slice by an element of the centre of \(SU(3)\). The Polyakov loop is not invariant, so \(\langle L\rangle\) is an order parameter of deconfinement. Dynamical quarks break this symmetry explicitly, so in full QCD \(L\) is only a probe, not an order parameter.
free energy vs energy	at \(T=0\) one talks about the static energy \(V_S(r)\) (an eigenvalue of the Hamiltonian). At \(T>0\) what the correlator gives is a free energy \(F = E - TS\): it contains the entropy of the medium rearranging itself around the pair.
colour singlet / octet	two static charges in the fundamental and antifundamental form \(3\otimes\bar 3 = 1\oplus 8\). The singlet channel is attractive (\(V_s = -C_F\alpha_s/r\)), the octet channel is repulsive.
screening	in a plasma a colour charge polarises the medium; the Coulomb tail \(1/r\) is replaced by a Yukawa tail \(e^{-m_D r}/r\). \(m_D\) is the Debye (electric screening) mass.
the three thermal scales	\(\pi T\) (hard, the lowest non-zero Matsubara frequency), \(gT\) (soft, electric / Debye), \(g^2T\) (ultrasoft, magnetic). Because \(g<1\) they are hierarchically separated at very high \(T\).
pNRQCD	potential non-relativistic QCD: the EFT you get after integrating out the scale \(1/r\). Valid when \(1/r\) is the largest scale, i.e. \(r \ll 1/T\). Its degrees of freedom are a singlet field and an octet field with potentials \(V_s, V_o\).
EQCD	electrostatic QCD: the 3-dimensional effective theory you get after integrating out the hard scale \(\pi T\). Its fields are the 3d gauge field and an adjoint scalar \(A_0\) with mass \(m_D\). Valid for \(r \sim 1/(gT)\).
Linde problem	at the magnetic scale \(g^2T\) the perturbative series stops converging; nothing there can be computed by expansion in \(g\).
HISQ	Highly Improved Staggered Quark action: a lattice fermion discretisation whose cutoff effects start at \(\mathcal{O}(a^2)\) and contain only even powers of \(a\). This is why continuum fits below use \(1/N_\tau^2, 1/N_\tau^4\).
\(N_\sigma^3 \times N_\tau\)	the lattice: \(N_\tau\) sites in the time direction, so \(T = 1/(a N_\tau)\). Taking \(N_\tau\to\infty\) at fixed \(T\) is the continuum limit.
\(\beta = 10/g_0^2\)	the lattice bare gauge coupling parameter. Large \(\beta\) = weak bare coupling = fine lattice.

\[ C_P(r,T) = \exp\left(-\frac{F_{Q\bar Q}(r,T)}{T}\right) \]

Words. \(C_P\) = renormalised correlator of two Polyakov loops. \(F_{Q\bar Q}\) = free energy of a static quark-antiquark pair separated by \(r\) at temperature \(T\).

Derivation. The thermal partition function is \(Z = \mathrm{Tr}\,e^{-H/T}\) and the free energy is \(F = -T\ln Z\). Inserting a static quark at \(\mathbf{x}\) and a static antiquark at \(\mathbf{x}+\mathbf{r}\) means restricting the trace to the sector of the Hilbert space containing that pair. In the Euclidean path integral that sector is produced exactly by inserting the two Polyakov loops, because a Polyakov loop is the propagator of an infinitely heavy quark for one full period \(1/T\). Hence

\[ \langle P(\mathbf{x}) P^\dagger(\mathbf{x}+\mathbf{r})\rangle = \frac{Z_{Q\bar Q}(r,T)}{Z(T)} = \frac{e^{-F_{Q\bar Q\,\text{total}}/T}}{e^{-F_0/T}} = e^{-\left(F_{Q\bar Q\,\text{total}}-F_0\right)/T} \]

and \(F_{Q\bar Q} \equiv F_{Q\bar Q\,\text{total}} - F_0\) is by definition the extra free energy caused by the pair. Taking logs gives Eq. (1). Note this is an equation of definition once \(C_P\) is renormalised; the renormalisation constant is fixed by demanding \(F_{Q\bar Q}(r\to\infty) = 2F_Q\).

\[ P(N_\tau,\mathbf{x}) = \frac{1}{3}\,\mathrm{Tr}\,W(N_\tau,\mathbf{x}), \qquad W(N_\tau,\mathbf{x}) = \prod_{\tau/a=1}^{N_\tau} U_0(\tau,\mathbf{x}) \]

Words. \(U_0(\tau,\mathbf{x})\in SU(3)\) is the temporal link variable, the lattice version of \(\exp\left(iag A_0(\tau,\mathbf{x})\right)\). \(W\) is the temporal Wilson line, \(N_\tau\) links long, i.e. exactly once around the periodic time direction.

Derivation. The continuum Polyakov loop is

\[ W(\mathbf{x}) = \mathcal{P}\exp\left(ig\int_0^{1/T} d\tau\, A_0(\tau,\mathbf{x})\right). \]

Discretise: split \([0,1/T]\) into \(N_\tau\) steps of length \(a\), so \(1/T = aN_\tau\). Each infinitesimal factor \(\exp(iagA_0)\) is one link \(U_0\), and path ordering becomes ordered matrix multiplication. That is Eq. (3). The \(1/3 = 1/N_c\) in Eq. (2) is a normalisation convention: with it, \(P = 1\) for the trivial configuration \(U_0 = \mathbb{1}\). This convention is the reason the constant \(T\ln 9 = T\ln N_c^2\) appears later.

\[ L_{\text{bare}}(\beta,N_\tau) = \langle P(N_\tau,\mathbf{x})\rangle \]

Expectation value, additionally averaged over all spatial sites \(\mathbf{x}\) (translation invariance, done to reduce noise). "bare" = not yet renormalised; it still carries a divergence.

\[ C_P^{\text{bare}}(\beta,N_\tau,r) = \langle P(N_\tau,\mathbf{x}) P^\dagger(N_\tau,\mathbf{x}+\mathbf{r})\rangle \]

Words. \(P^\dagger\) is the antiquark (\(\bar 3\) transforms with the conjugate matrix). The average is over gauge configurations and over \(\mathbf{x}\).

\[ C_S^{\text{bare}}(\beta,N_\tau,r) = \frac{1}{3}\left\langle \mathrm{Tr}\left[W(N_\tau,\mathbf{x})W^\dagger(N_\tau,\mathbf{x}+\mathbf{r})\right]\right\rangle \qquad \text{(Coulomb gauge)} \]

Words. Here the two Wilson lines are inside one trace: the colour indices of quark and antiquark are contracted with each other, which is precisely the colour-singlet projection. In Eq. (5) they were traced separately, which is a colour average over singlet and octet.

Derivation of why this is "the singlet". Write the two static-quark states with colour indices \(a,b\). The singlet projector is \(\delta^{ab}/\sqrt{N_c}\). The standard McLerran-Svetitsky decomposition is

\[ e^{-F_1/T} = \frac{1}{N_c}\left\langle \mathrm{Tr}\left[W(\mathbf{x})W^\dagger(\mathbf{y})\right]\right\rangle, \qquad e^{-F_8/T} = \frac{1}{N_c^2-1}\left\langle \mathrm{Tr}W\,\mathrm{Tr}W^\dagger\right\rangle - \frac{1}{N_c(N_c^2-1)}\left\langle \mathrm{Tr}\left[WW^\dagger\right]\right\rangle . \]

Check the consistency with Eq. (5). Since \(3\otimes\bar3 = 1\oplus 8\) has \(1 + 8 = 9 = N_c^2\) states,

\[ \frac{1}{N_c^2}e^{-F_1/T} + \frac{N_c^2-1}{N_c^2}e^{-F_8/T} = \frac{1}{N_c^3}\langle \mathrm{Tr}WW^\dagger\rangle + \frac{1}{N_c^2}\langle \mathrm{Tr}W\,\mathrm{Tr}W^\dagger\rangle - \frac{1}{N_c^3}\langle \mathrm{Tr}WW^\dagger\rangle = \langle P P^\dagger\rangle . \]

The \(\mathrm{Tr}WW^\dagger\) pieces cancel exactly and one recovers the Polyakov loop correlator. So Eq. (6) is the singlet piece of Eq. (5). This algebra is the seed of Eq. (26) below.

Why Coulomb gauge? \(\mathrm{Tr}[W(\mathbf{x})W^\dagger(\mathbf{y})]\) with \(\mathbf{x}\ne\mathbf{y}\) is not gauge invariant: strictly one would need a spatial Wilson line joining the two ends. Fixing to Coulomb gauge is the standard substitute; it is the gauge in which the object stays UV finite (in covariant gauges the normalised correlator is UV divergent) and in which weak-coupling results exist.

\[ W_S^{\text{bare}}(\beta,N_\tau,r) = \frac{1}{3}\left\langle \mathrm{Tr}\left[W(N_\tau,\mathbf{x})S(N_\tau,\mathbf{x},\mathbf{r})W^\dagger(N_\tau,\mathbf{x}+\mathbf{r})S^\dagger(0,\mathbf{x},\mathbf{r})\right]\right\rangle \]

Words. \(S(\tau,\mathbf{x},\mathbf{r})\) is a spatial Wilson line from \((\tau,\mathbf{x})\) to \((\tau,\mathbf{x}+\mathbf{r})\). Putting the two spatial lines in closes the contour, so this object is gauge invariant without gauge fixing. "Cyclic" = the temporal sides wrap the whole period, so the contour is smooth (no corners) in the continuum.

Price paid. The closed contour brings extra divergences: self-energy of the spatial lines, intersection divergences where spatial and temporal lines meet, and, on the lattice, cusp divergences for off-axis \(\mathbf{r}\) (e.g. \(r/a=\sqrt2\)). Hence the paper uses only on-axis separations and applies HYP smearing to the spatial links to shrink the self-energy divergence.

The central fact: a Wilson line of length \(\ell\) has a linearly divergent self-energy, \(\propto e^{-\delta m\,\ell}\) with \(\delta m \sim c/a\). A Polyakov loop has \(\ell = 1/T\), so

\[ L_{\text{bare}} = e^{-C_Q/T}\,L, \qquad C_Q \xrightarrow[a\to0]{} \infty . \]

The correlators (5) and (6) contain two such lines, hence exactly the square of that factor. Every ratio in which the divergence appears twice upstairs and twice downstairs is therefore finite.

\[ C_P^{\text{sub}}(\beta,N_\tau,r) = \frac{C_P^{\text{bare}}(\beta,N_\tau,r)}{\left[L_{\text{bare}}(\beta,N_\tau)\right]^2} = \frac{C_P(\beta,N_\tau,r)}{\left[L(\beta,N_\tau)\right]^2} \]

Derivation. Substituting \(C_P^{\text{bare}} = e^{-2C_Q/T}C_P\) and \(L_{\text{bare}}^2 = e^{-2C_Q/T}L^2\), the divergent factors cancel identically, giving the second equality. So \(C_P^{\text{sub}}\) has a continuum limit and is scheme independent.

Words. "sub" = subtracted (in the free energy it will be a subtraction, since dividing correlators = subtracting free energies). \(C_P^{\text{sub}}\) is the connected part of the correlator: at \(r\to\infty\) clustering gives \(C_P^{\text{bare}}\to L_{\text{bare}}^2\), so \(C_P^{\text{sub}}\to 1\). All information about the mutual interaction sits in the deviation of \(C_P^{\text{sub}}\) from 1; the disconnected part \(L^2\) describes two quarks interacting only with the medium.

\[ C_S^{\text{sub}}(\beta,N_\tau,r) = \frac{C_S^{\text{bare}}(\beta,N_\tau,r)}{\left[L_{\text{bare}}(\beta,N_\tau)\right]^2} = \frac{C_S(\beta,N_\tau,r)}{\left[L(\beta,N_\tau)\right]^2} \]

Same argument: \(C_S^{\text{bare}}\) also contains exactly two temporal lines and (in Coulomb gauge) no spatial line, so its divergence is again \(e^{-2C_Q/T}\). In a covariant gauge this statement fails.

\[ F^{\text{sub}}_{Q\bar Q}(r,T,a) = -T\ln C_P^{\text{sub}}(\beta,N_\tau,r), \qquad F^{\text{sub}}_S(r,T,a) = -T\ln C_S^{\text{sub}}(\beta,N_\tau,r) \]

Definition, by analogy with Eq. (1). The trade \((\beta,N_\tau)\to(T,a)\) uses \(T = 1/(N_\tau a)\) and the known \(a(\beta)\).

\[ L(T) = \exp\left(-\frac{F_Q(T)}{T}\right) \]

Derivation. Same argument as Eq. (1) but with a single Polyakov loop: \(\langle P\rangle = Z_Q/Z = e^{-F_Q/T}\), where \(F_Q\) is the free energy cost of putting one static quark into the medium.

\[ F_{Q\bar Q} = F^{\text{sub}}_{Q\bar Q} + 2F_Q, \qquad F_S = F^{\text{sub}}_S + 2F_Q \]

Derivation. Combine (10), (8), (1), (12):

\[ F^{\text{sub}}_{Q\bar Q} = -T\ln\frac{C_P}{L^2} = -T\ln C_P + 2T\ln L = F_{Q\bar Q} + 2\left(-F_Q\right), \]

using \(T\ln L = -F_Q\). Rearranging gives (13). Eq. (14) is the same statement for the singlet channel (there it is a definition of \(F_S\), since no independent equation fixes the singlet normalisation).

Reading. \(-2T\ln L = 2F_Q\) is the free energy of two static charges at infinite separation, each dressed by the medium but not talking to each other. So the whole renormalisation problem of the pair reduces to the renormalisation of one Polyakov loop.

*The \(T\ln 9\).* Because of the \(1/N_c\) in Eq. (2), \(C_P\) carries an extra \(1/N_c^2 = 1/9\) relative to the usual convention. At short distance \(C_P\to \tfrac19 e^{-V_s/T}\), so \(F_{Q\bar Q}\to V_s + T\ln 9\). That is why figures show \(F_{Q\bar Q}-T\ln 9\) when comparing to \(V_S\).

\[ \alpha_{Q\bar Q}(r) = \frac{1}{C_F}\,r^2\,\frac{\partial V_S(r)}{\partial r}, \qquad C_F = \frac{N_c^2-1}{2N_c} = \frac43 \]

Words. \(C_F\) = quadratic Casimir of the fundamental representation of \(SU(3)\); it is the colour factor of one-gluon exchange between a \(3\) and a \(\bar3\) in the singlet channel. \(V_S(r)\) = static energy at \(T=0\). "Effective coupling" = a coupling defined from a physical quantity, so it is scheme independent (unlike \(\alpha_s^{\overline{\rm MS}}\)).

Derivation. At tree level \(V_S(r) = -C_F\alpha_s/r\), hence \(\partial V_S/\partial r = C_F\alpha_s/r^2\) and

\[ \frac{1}{C_F}r^2\frac{\partial V_S}{\partial r} = \alpha_s . \]

So Eq. (15) is constructed so that it equals \(\alpha_s\) at leading order and, beyond leading order, defines a running coupling in the "force scheme" with \(\mu\sim 1/r\).

Use in the paper. The same formula is applied with \(V_S \to F_S\) and \(V_S\to F_{Q\bar Q}\). At short \(r\) the coupling runs like the vacuum one; it reaches a maximum at \(r_{\max}(T)\approx 0.4/T\) and then falls -- the fall is the onset of screening. \(\alpha_{Q\bar Q}(r_{\max},T)<0.5\) only for \(T\gtrsim320\) MeV, which is where a weak-coupling treatment can be hoped for.

Setting. \(rT \ll 1\), i.e. \(1/r \gg T\). Within that, choose the sub-hierarchy \(1/r \gg T \gg \alpha_s/r\) and \(T\gg m_D\). Then the pair is essentially a vacuum object and the medium is a small perturbation. The paper compares \(V_S - F_S\) (the difference removes the hard-to-compute vacuum potential and cancels lattice artefacts and the additive renormalisation).

\[ V_S - F_S = \left(\Delta F_g + \Delta F_f + \Delta F_S\right)T + \mathcal{O}\left(g^4(rT)^5,\,g^6\right) \]

Words. The \(\Delta F\)'s are dimensionless; the overall \(T\) carries the dimension, because every thermal correction must vanish as \(T\to0\). Three sources: \(\Delta F_g\) from non-static gluons (Matsubara modes \(\omega_n = 2\pi nT\), \(n\ne0\)), \(\Delta F_f\) from \(N_f\) massless quarks in the loop, \(\Delta F_S\) from the soft/Debye scale \(m_D\).

\[ \Delta F_g = \alpha_s^2 N_c C_F\left[-\frac{\pi}{9}(rT) + \frac43\zeta(3)(rT)^2 - \frac{22\pi^3}{675}(rT)^3\right] \] \[ \Delta F_f = \alpha_s^2 N_f C_F\left[\zeta(3)(rT)^2 - \frac{7\pi^3}{270}(rT)^3\right] \]

Structure (this is what to understand, the coefficients are a two-loop computation from Ref. [31]).

Both are \(\mathcal{O}(\alpha_s^2)=\mathcal{O}(g^4)\): they come from two powers of the coupling, i.e. the dipole \(\times\) dipole insertion of the thermal gluon/quark loop on the static pair.
The colour factors are transparent: \(C_F\) is the colour charge of the static pair, \(N_c\) the gluon-loop colour factor, \(N_f\) counts quark loops.
The expansion parameter is \(rT\): this is the multipole expansion of pNRQCD. The pair of size \(r\) probes a thermal field of wavelength \(1/T\); the leading coupling is dipole, \(\sim r\), and each further power of \(r\) costs one more power of \(rT\).
\(\zeta(3)\) and \(\pi^3\) appear because the Matsubara sums \(\sum_{n\neq0} n^{-k}\) produce \(\zeta\) values, and \(\zeta(2)=\pi^2/6\), \(\zeta(4)=\pi^4/90\) etc.
The quark contribution has no linear term in \(rT\). Fermions obey antiperiodic boundary conditions in \(\tau\), so their Matsubara frequencies are the half-integer ones \(\omega_n = (2n+1)\pi T\); the sum structure that produces the \((rT)^1\) piece for gluons has no fermionic counterpart.

\[ \Delta F_S = -\alpha_s^{5/2}\,\frac{C_F}{6}\left[4\pi\,\frac{N_c+\frac{N_f}{2}}{3}\right]^{3/2}(rT)^2 \]

This one can be derived cleanly. The soft scale contributes through the screened Coulomb potential. Write the resummed singlet free energy at short distance,

\[ F_S = -C_F\alpha_s\frac{e^{-m_D r}}{r} - C_F\alpha_s m_D . \]

(The constant \(-C_F\alpha_s m_D\) is the self-energy of each charge in the plasma; it is what makes \(F_S\to\) finite at \(r\to0\).) Expand for \(m_D r\ll1\):

\[ F_S = -\frac{C_F\alpha_s}{r}\left(1 - m_Dr + \frac{m_D^2r^2}{2} - \frac{m_D^3r^3}{6}+\dots\right) - C_F\alpha_s m_D = -\frac{C_F\alpha_s}{r} - \frac{C_F\alpha_s m_D^2 r}{2} + \frac{C_F\alpha_s m_D^3 r^2}{6}+\dots \]

The \(\mathcal{O}(m_D)\) terms cancel against the constant. Since \(V_S = -C_F\alpha_s/r\) at this order,

\[ V_S - F_S = +\frac{C_F}{2}\alpha_s m_D^2 r - \frac{C_F}{6}\alpha_s m_D^3 r^2 + \dots \]

*Which of these belongs in \(\Delta F_S\)?* Only the odd powers of \(m_D\). Since \(m_D\propto g T\), an even power \(m_D^{2k}\) is an integer power of \(\alpha_s\) and is therefore already contained in the strict fixed-order (hard) calculation, i.e. in \(\Delta F_g+\Delta F_f\) -- indeed \(\alpha_s m_D^2 r \sim \alpha_s^2 (rT) T\), which is exactly the \((rT)^1\) structure of Eq. (17). An odd power \(m_D^{2k+1}\) is non-analytic in \(\alpha_s\) (a half-integer power) and can only come from the resummed soft region; it is genuinely new. The first such term is the \(m_D^3 r^2\) one:

\[ V_S - F_S \supset -\frac{C_F}{6}\alpha_s m_D^3 r^2 , \]

which is why Eq. (19) contains a \((rT)^2\) term and nothing else.
Now insert the leading-order Debye mass (Eq. 20 below) in the form

\[ \frac{m_D}{T} = \sqrt{4\pi\alpha_s\,\frac{N_c+\frac{N_f}{2}}{3}} \quad\Longrightarrow\quad \left(\frac{m_D}{T}\right)^3 = \alpha_s^{3/2}\left[4\pi\frac{N_c+\frac{N_f}{2}}{3}\right]^{3/2}, \]

so that

\[ V_S-F_S \supset -\frac{C_F}{6}\alpha_s T^3 \left(\frac{m_D}{T}\right)^3 r^2 = -\alpha_s^{5/2}\frac{C_F}{6}\left[4\pi\frac{N_c+\frac{N_f}{2}}{3}\right]^{3/2}(rT)^2\,T , \]

which is exactly \(\Delta F_S\,T\) with \(\Delta F_S\) as written. The half-integer power \(\alpha_s^{5/2}\), i.e. order \(g^5\), is the signature of a soft, resummed contribution: odd powers of \(g\) can only come from the \(m_D\propto gT\) scale.

Physics read-off.

\(V_S-F_S \to 0\) as \(r\to0\): at zero separation the pair is colour neutral and point-like, the medium cannot resolve it. The lattice data instead approach a non-zero constant in units of \(T\); this is attributed to an \(r\)-independent \(\mathcal{O}(g^6)T\) term from the pNRQCD-to-NRQCD matching that is not in Eq. (16). That is why Fig. 7 shifts the perturbative curve by a constant.
The \((rT)^2\) terms of (17)+(18) are positive, that of (19) is negative: a partial cancellation. This is exactly why the measured medium effect on \(F_S\) is so small out to \(rT\approx0.3\).
\(\Delta F_S\) is parametrically smaller (\(g^5\) vs \(g^4\)) but numerically dominant once \(rT\) grows -- the signal that the hierarchy \(1/r\gg T\gg m_D\) has broken down and one must resum, i.e. go to EQCD.

\[ m_D\big|_{\rm LO}(\mu) = g(\mu)\,T\sqrt{\frac{2N_c+N_f}{6}} \]

Words. \(m_D\) is the mass of the adjoint scalar \(A_0\) in EQCD, equivalently the inverse screening length of chromoelectric fields.

Derivation sketch. Compute the one-loop gluon self-energy \(\Pi_{00}(k_0=0,\mathbf{k}\to0)\) at finite \(T\). Gluons contribute \(g^2T^2N_c/3\) and each massless Dirac flavour \(g^2T^2/6\), so

\[ m_D^2 = \Pi_{00}(0,0) = g^2T^2\left(\frac{N_c}{3}+\frac{N_f}{6}\right) = g^2T^2\,\frac{2N_c+N_f}{6}. \]

The \(T^2\) is fixed by dimensions (the only scale is \(T\)), and \(m_D\propto gT\) places it at the soft scale. Equivalently \(m_D^2 = 4\pi\alpha_s T^2(N_c+N_f/2)/3\), the form used in Eq. (19).

\[ F^{\text{sub}}_S\big|_{\rm NLO}(r,T) = -C_F\frac{\alpha_s(\mu)}{r}e^{-rm_D(\mu)} \Big(1+\alpha_s(\mu)N_c rT\left[2-\ln\left(2rm_D(\mu)\right)-\gamma_E + e^{+2rm_D(\mu)}E_1\left(2rm_D(\mu)\right)\right]\Big) \]

Words. This is Eq. (3.22) of Ref. [42], with the two-gluon-exchange piece \(\propto g^4/r^2\) dropped (in the regime \(r\sim1/m_D\) that piece is \(\mathcal{O}(g^6)\), beyond NLO).

How to read it.

The prefactor \(-C_F\alpha_s e^{-m_Dr}/r\) is the Debye-screened Coulomb potential: exchange of one \(A_0\) field of mass \(m_D\) in three dimensions gives the Yukawa propagator, and \(\int\!\frac{d^3k}{(2\pi)^3}\frac{e^{i\mathbf{k}\cdot\mathbf{r}}}{\mathbf{k}^2+m_D^2}=\frac{e^{-m_Dr}}{4\pi r}\).
The bracket is the NLO correction, one more power of \(\alpha_s\) times \(N_c\) (a gluon self-interaction) times \(rT\) (the 3d loop expansion parameter of EQCD is \(g^2T/m_D\cdot\)..., which for a line of length \(r\) produces \(\alpha_s N_c rT\)).
The logarithm \(\ln(2rm_D)\) and \(\gamma_E\) come from the 3d one-loop integral; they signal that the NLO piece is sensitive to both scales \(1/r\) and \(m_D\).
\(E_1\) makes the expression well behaved at both ends: the combination \(e^{+2rm_D}E_1(2rm_D)\) tends to \(-\ln(2rm_D)-\gamma_E+\dots\) for small argument (cancelling part of the explicit log) and to \(1/(2rm_D)\) for large argument.

\[ E_1(z) = \int_z^\infty \frac{dt}{t}e^{-t} \]

The exponential integral. Small-\(z\) behaviour \(E_1(z) = -\gamma_E - \ln z + z - \dots\); large-\(z\) behaviour \(E_1(z)\simeq e^{-z}/z\).

\[ \delta F^{\text{sub}}_S(r,T) = -\frac{C_F\alpha_s(\mu)}{r}e^{-rm_D(\mu)}\left(1-\frac{rm_D(\mu)}{2}\right)\delta Z_1(\mu) \] \[ \delta Z_1(\mu) = \alpha_s(\mu)\left[\frac{11N_c}{3}+\frac23\left(1-4\ln2\right)+2\beta_0\left(\gamma_E+\ln\frac{\mu}{4\pi T}\right)\right], \qquad \beta_0 = \frac{11N_c}{3}-\frac{2N_f}{3} \]

Words. \(\beta_0\) is the one-loop beta-function coefficient, \(\mu\,d\alpha_s/d\mu = -\beta_0\alpha_s^2/(2\pi)+\dots\). \(\delta Z_1\) is a matching / normalisation constant of EQCD.

Where it comes from. When the hard scale \(\pi T\) is integrated out to build EQCD, the 3d coupling and the \(A_0\) field pick up matching coefficients. The piece \(2\beta_0\left(\gamma_E+\ln\frac{\mu}{4\pi T}\right)\) is precisely what converts \(\alpha_s\) evaluated at a generic \(\mu\) into \(\alpha_s\) evaluated at the natural thermal scale: the combination \(\mu\,e^{\gamma_E}/(4\pi T)\) inside a log times \(\beta_0\) is the standard one-loop running. The remaining \(\frac{11N_c}{3}+\frac23(1-4\ln2)\) is a pure constant from the hard one-loop matching (\(\ln 2\)'s are typical of Bose/Fermi Matsubara sums). The structure \(\left(1-\frac{rm_D}{2}\right)\) arises from differentiating the screened potential with respect to the coupling that also sits inside \(m_D\propto g\):

\[ \frac{\partial}{\partial \ln g}\left(-C_F\frac{\alpha_s}{r}e^{-m_Dr}\right)\ \propto\ -C_F\frac{\alpha_s}{r}e^{-m_Dr}\left(2 - m_D r\right)\ \propto\ \left(1-\frac{rm_D}{2}\right). \]

So Eq. (23) is "the effect of letting the coupling run", formally higher order but numerically important because \(\beta_0\) is large.

\[ m_D^2\big|_{\rm NLO}(\mu) = m_D^2\big|_{\rm LO}(\mu)\left(1+\frac{\alpha_s(\mu)}{4\pi}\left[2\beta_0\left(\gamma_E+\ln\frac{\mu}{4\pi T}\right)+\frac{5N_c}{3}+\frac{2N_f}{3}\left(1-4\ln2\right)\right]\right)-C_FN_f\alpha_s^2(\mu)T^2 \]

Reading. Three ingredients:

the \(\beta_0\) log, again just the running of \(g^2\) inside \(m_D^2 = g^2T^2(2N_c+N_f)/6\) from \(\mu\) to \(\sim 4\pi T e^{-\gamma_E}\);
genuine two-loop constants \(\frac{5N_c}{3}\) and \(\frac{2N_f}{3}(1-4\ln2)\) from the hard matching;
a separate, non-multiplicative term \(-C_FN_f\alpha_s^2T^2\): a quark-loop contribution that is not proportional to the LO structure, hence written outside.

Note \(m_D^2|_{\rm NLO}\) is well defined (gauge invariant) precisely because it is a parameter of EQCD, fixed by matching, rather than a pole of a gauge-dependent propagator. All numerical work in the paper uses this NLO \(m_D\).

\[ C_P(r,T)\equiv e^{-F_{Q\bar Q}/T} = \frac{1}{N_c^2}e^{-f_s/T} + \frac{N_c^2-1}{N_c^2}e^{-f_o/T} \]

Words. \(f_s\), \(f_o\) are gauge-invariant singlet and octet free energies defined through pNRQCD correlators (small letters = perturbative/EFT objects; capital letters \(F_S\), \(F_O\) = lattice objects; they need not coincide).

Derivation. This is the algebra already done under Eq. (6), read backwards. \(3\otimes\bar3 = 1\oplus8\) has \(N_c^2\) states in total, \(1\) singlet and \(N_c^2-1\) octet. The Polyakov loop correlator is the colour-averaged (unpolarised in colour) object, so the two channels enter with their multiplicity weights \(1/N_c^2\) and \((N_c^2-1)/N_c^2\). Each channel Boltzmann-weights with its own free energy.

Immediate consequence. \(-T\ln\) of a sum of two exponentials is not the sum of the two free energies. The Polyakov loop correlator is therefore never simply "the potential", and the octet piece -- repulsive, and weighted 8/9 -- dominates as soon as it is not suppressed.

\[ e^{-F_{Q\bar Q}/T} = \frac{1}{N_c^2}e^{-V_s/T} + \frac{N_c^2-1}{N_c^2}L_A\,e^{-V_o/T}, \qquad V_s = -\frac{C_F\alpha_s}{r}+\dots,\quad V_o = -\frac{V_s}{N_c^2-1}+\dots \]

Words. \(V_s,V_o\) are the \(T=0\) singlet and octet potentials. \(L_A\) is the adjoint Polyakov loop expectation value.

Derivation / logic. At \(r\ll1/T\) the pair is a vacuum object, so \(f_s\to V_s\), i.e. the singlet term is purely a vacuum quantity. In the octet channel the pair carries a net adjoint colour charge which must be screened by the medium; the cost of that is exactly the free energy of an adjoint static source, i.e. the factor \(L_A = e^{-F_A/T}\). This is why \(L_A\) multiplies the octet term. The relation \(V_o = -V_s/(N_c^2-1)\) is Casimir scaling at one gluon exchange: the colour factor of one-gluon exchange is \(-C_F\) in the singlet and \(+1/(2N_c) = C_F/(N_c^2-1)\) in the octet, i.e. the octet is repulsive and weaker by \(-(N_c^2-1)\). Eq. (27) is valid up to \(\mathcal{O}(\alpha_s^3)\).

Physical read-off. At low \(T\), \(L_A\) is tiny, so the octet is suppressed and \(F_{Q\bar Q}\approx V_s + T\ln N_c^2\). At high \(T\), \(L_A\sim1\) and the repulsive octet, with weight \(8/9\), takes over -- which is why \(F_{Q\bar Q}\) departs from \(V_S\) much earlier than \(F_S\) does.

\[ C_o = \frac{N_c^2-1}{N_c^2}e^{-f_o/T} \]

Just the second term of Eq. (26), given a name. On the lattice the proxy is \(C_O \equiv C_P - e^{-V_S/T}/N_c^2\), i.e. "measured correlator minus the vacuum singlet piece", and \(F_O = -T\ln\!\big[N_c^2 C_O/(N_c^2-1)\big]\) is the corresponding octet free energy.

\[ L_A = L^{C_A/C_F}\cdot\left(1-\delta_8\right)^{C_A/C_F}, \qquad C_F = \frac{N_c^2-1}{2N_c},\ \ C_A = N_c \]

Words. \(C_A\) = Casimir of the adjoint representation. \(\delta_8\) = a measured parameter quantifying the violation of Casimir scaling; \(\delta_8=0\) means exact Casimir scaling.

Derivation. In perturbation theory the free energy of a static source in representation \(R\) is proportional to its Casimir, \(F_R = (C_R/C_F)F_Q + \dots\). Hence

\[ L_R = e^{-F_R/T} = \left(e^{-F_Q/T}\right)^{C_R/C_F} = L^{C_R/C_F}. \]

For the adjoint \(C_A/C_F = 2N_c^2/(N_c^2-1) = 9/4\) for \(N_c=3\). The factor \((1-\delta_8)^{C_A/C_F}\) is the empirical correction; lattice studies find \(\delta_8\) small above \(T_c\) and consistent with zero for \(T>300\) MeV.

\[ e^{-F^{\text{sub}}_{Q\bar Q}/T} = \frac{1}{N_c^2}\exp\left[-\frac{V_S-2F_Q}{T}\right] +\frac{N_c^2-1}{N_c^2}\exp\left[\frac{V_S-2F_Q}{(N_c^2-1)T}\right]\cdot\left(1-\delta_8\right)^{C_A/C_F}\cdot e^{-\delta V_o/T} \] \[ \delta V_o = V_o + \frac{V_s}{N_c^2-1} \]

Derivation (do this one carefully; it is the paper's main quantitative test). Start from Eq. (27) and multiply both sides by \(e^{2F_Q/T}\), using \(F^{\text{sub}}_{Q\bar Q} = F_{Q\bar Q}-2F_Q\) from Eq. (13):

Singlet term: \(\frac{1}{N_c^2}e^{-V_s/T}e^{2F_Q/T} = \frac{1}{N_c^2}e^{-(V_s-2F_Q)/T}\).
Octet term. Write \(V_o = -\frac{V_s}{N_c^2-1}+\delta V_o\) (this defines \(\delta V_o\) as the Casimir-scaling-violating remainder), so \(e^{-V_o/T}=e^{+V_s/((N_c^2-1)T)}e^{-\delta V_o/T}\).
Use Eq. (29) with \(L = e^{-F_Q/T}\): \(L_A = e^{-\frac{C_A}{C_F}\frac{F_Q}{T}}(1-\delta_8)^{C_A/C_F}\), and \(\frac{C_A}{C_F} = \frac{2N_c^2}{N_c^2-1}\).
Collect the \(F_Q\) exponents in the octet term: \[ \frac{2F_Q}{T}-\frac{C_A}{C_F}\frac{F_Q}{T} = \frac{2F_Q}{T}\left(1-\frac{N_c^2}{N_c^2-1}\right) = -\frac{2F_Q}{(N_c^2-1)T}. \]
Combine with step 2: \[ e^{\frac{V_s}{(N_c^2-1)T}}\,e^{-\frac{2F_Q}{(N_c^2-1)T}} = \exp\left[\frac{V_s-2F_Q}{(N_c^2-1)T}\right]. \]

Putting 1-5 together gives Eq. (30) exactly, with the lattice \(T=0\) static energy \(V_S\) used as a proxy for \(V_s\).

Why it matters. Every ingredient on the right-hand side except \(\delta V_o\) is a lattice number (\(V_S\) from \(T=0\), \(F_Q\) from the renormalised Polyakov loop, \(\delta_8\) measured). So Eq. (30) predicts \(F^{\text{sub}}_{Q\bar Q}\) with essentially one perturbative input, and the comparison with the direct measurement (Fig. 10) tests pNRQCD. It works up to \(rT\approx0.3\).

\[ \delta V_o = \frac{N_c\,\alpha_s^3}{8\,r}\left(\frac{\pi^2}{4}-3\right)+\mathcal{O}(\alpha_s^4) \]

Words. At one gluon exchange the octet potential is exactly \(-V_s/(N_c^2-1)\). The first deviation appears at three loops in the colour-Coulomb sense, order \(\alpha_s^3\), from the non-abelian three-gluon vertex; \(N_c\) in the numerator is its colour factor and \(\pi^2/4-3\approx-0.533\) its (negative) coefficient, so \(\delta V_o<0\). It still falls like \(1/r\), so it is a genuine potential. Numerically it matters at high \(T\): without it the reconstruction of Eq. (30) only works up to \(rT<0.2\) instead of \(0.3\).

Eq. (32) : leading order \(F^{\text{sub}}_{Q\bar Q}\) for \(\alpha_s/r\ll T\)

\[ F^{\text{sub}}_{Q\bar Q} = -\frac{N_c^2-1}{8N_c^2}\frac{\alpha_s^2}{r^2T} \]

Full derivation. Take Eq. (26)/(27) with \(L_A\to1\) and expand both exponentials, since \(V_{s,o}/T\sim\alpha_s/(rT)\ll1\):

\[ C_P^{\text{sub}} \simeq \frac{1}{N_c^2}\left(1-\frac{V_s}{T}+\frac{V_s^2}{2T^2}\right)+\frac{N_c^2-1}{N_c^2}\left(1-\frac{V_o}{T}+\frac{V_o^2}{2T^2}\right). \]

\(\mathcal{O}(1)\): \(\frac{1}{N_c^2}+\frac{N_c^2-1}{N_c^2} = 1\). Good, \(C_P^{\text{sub}}\to1\).
\(\mathcal{O}(1/T)\), with \(V_s = -\frac{C_F\alpha_s}{r}\) and \(V_o = +\frac{C_F\alpha_s}{(N_c^2-1)r}\): \[ -\frac{1}{T}\left[\frac{1}{N_c^2}\left(-\frac{C_F\alpha_s}{r}\right)+\frac{N_c^2-1}{N_c^2}\cdot\frac{C_F\alpha_s}{(N_c^2-1)r}\right] = -\frac{1}{T}\left[-\frac{C_F\alpha_s}{N_c^2 r}+\frac{C_F\alpha_s}{N_c^2 r}\right]=0 . \]

*This exact cancellation of the \(1/r\) term between singlet and octet is the central point*: the leading Coulomb behaviour disappears from the Polyakov loop correlator.

\(\mathcal{O}(1/T^2)\): \[ \frac{1}{2T^2}\left[\frac{1}{N_c^2}\frac{C_F^2\alpha_s^2}{r^2}+\frac{N_c^2-1}{N_c^2}\frac{C_F^2\alpha_s^2}{(N_c^2-1)^2r^2}\right] =\frac{C_F^2\alpha_s^2}{2N_c^2r^2T^2}\left[1+\frac{1}{N_c^2-1}\right] =\frac{C_F^2\alpha_s^2}{2r^2T^2\left(N_c^2-1\right)} . \]

Using \(C_F^2 = \frac{(N_c^2-1)^2}{4N_c^2}\) this is \(\frac{N_c^2-1}{8N_c^2}\frac{\alpha_s^2}{r^2T^2}\). So

\[ C_P^{\text{sub}} = 1+\frac{N_c^2-1}{8N_c^2}\frac{\alpha_s^2}{r^2T^2}+\dots \ \Longrightarrow\ F^{\text{sub}}_{Q\bar Q}=-T\ln C_P^{\text{sub}} \simeq -\frac{N_c^2-1}{8N_c^2}\frac{\alpha_s^2}{r^2T}, \]

which is Eq. (32), valid up to \(\mathcal{O}(g^5)\). For \(N_c=3\) the prefactor is \(8/72 = 1/9\), whence the combination \(-9r^2TF^{\text{sub}}_{Q\bar Q}\) plotted in Figs. 11-13.

Read-off. \(F^{\text{sub}}_{Q\bar Q}\sim -1/r^2\), not \(-1/r\): the free energy of the pair is not Coulombic. Physically the surviving term is two-gluon exchange, the first process that is not cancelled between the channels. This is exactly the criterion used in Fig. 3: \(-r^2TF^{\text{sub}}_{Q\bar Q}\) constant \(\Rightarrow\) \(1/r^2\) regime (\(T\gg\alpha_s/r\)); \(-r^2TF^{\text{sub}}_{Q\bar Q}\) rising linearly \(\Rightarrow\) \(1/r\) regime (\(T\ll\alpha_s/r\), singlet-dominated).

\[ F^{\text{sub}}_{Q\bar Q} = -\frac{N_c^2-1}{8N_c^2}\frac{\alpha_s^2}{r^2T}\exp\left(-2m_Dr\right) \]

Derivation. Repeat the computation of Eq. (32), but with each of the two exchanged gluons being an electrostatic (\(A_0\)) gluon in the medium, so each propagator carries \(e^{-m_Dr}\) instead of being unscreened. Two exchanges \(\Rightarrow\) \(\left(e^{-m_Dr}\right)^2 = e^{-2m_Dr}\).

Consequences.

\(F^{\text{sub}}_{Q\bar Q}\) decays with twice the rate of \(F^{\text{sub}}_S\) (which decays as \(e^{-m_Dr}\)). Hence \(|F^{\text{sub}}_{Q\bar Q}|\ll|F^{\text{sub}}_S|\) at large \(r\), and the Polyakov loop correlator has a dreadful signal-to-noise ratio there. That is why the paper needs \(2\)-\(6\times10^4\) configurations and large aspect ratios \(N_\sigma/N_\tau=6\).
The measured screening masses satisfy \(m_{Q\bar Q}\approx2m_S\) at high \(T\), as Figs. 14-15 show, with \(m_S\approx(1.6\text{-}2.0)m_D\) and \(m_{Q\bar Q}\approx1.25\times2m_D\): the naive perturbative expectation is right up to a \(\mathcal{O}(1)\) non-perturbative enhancement coming from magnetic (\(g^2T\)) physics.

\[ \Delta f^{\text{bare}}_Q = \frac{f^{\text{bare}}_Q(m_1)-f^{\text{bare}}_Q(m_2)}{\sqrt{\left(\delta f^{\text{bare}}_Q(m_1)\right)^2+\left(\delta f^{\text{bare}}_Q(m_2)\right)^2}} \]

Words. This is a pull (a difference expressed in units of its own error). \(f_Q = F_Q/T\), \(\delta f\) = statistical error, \(m_1 = m_s/20\) and \(m_2 = m_s/5\) are two light sea-quark masses.

Derivation. Two independent measurements \(A\pm\delta A\) and \(B\pm\delta B\): the difference \(A-B\) has variance \(\delta A^2+\delta B^2\) (errors add in quadrature for independent samples). Dividing by that standard deviation makes the quantity \(\mathcal{N}(0,1)\)-distributed if the two agree. The paper finds \(|\Delta f^{\text{bare}}_Q|<1.5\), i.e. no significant quark-mass dependence, which justifies mixing the two mass sets.

\[ S_1 = -r\,F^{\text{sub}}, \qquad F\in\left\{F_S,\ F_{Q\bar Q}\right\} \]

Why this combination. If \(F^{\text{sub}} = -A\frac{e^{-mr}}{r}\), then \(S_1 = Ae^{-mr}\) is a pure exponential. Plotting \(\ln S_1\) against \(r\) therefore gives a straight line whose slope is minus the screening mass, with the Coulomb-like \(1/r\) prefactor divided out. This is the standard way to read a screening mass off a plot.

\[ S_Q = -\frac{\partial F_Q}{\partial T} \]

Standard thermodynamics: from \(F = E-TS\) and \(dF = -SdT\) at fixed volume. \(S_Q\) is the entropy the medium gains by accommodating one static colour charge. It is a sharper probe than \(F_Q\) (a derivative sharpens structure) and the paper shows it matches the NNLO weak-coupling result for \(T\gtrsim1.7\) GeV.

The problem. On a lattice, a separation of \(r/a = 1,\sqrt2,\sqrt3,2,\dots\) lattice steps does not see a rotationally symmetric world. The measured \(V_S(r)\) therefore wiggles at small \(r/a\): this is breaking of rotational symmetry, a pure discretisation artefact, and it is much larger than the statistical errors.

\[ V_S^{\text{lat,free}}(r) = \int\frac{d^3k}{(2\pi)^3}D_{00}(k_0=0,\mathbf{k})\,e^{i\mathbf{k}\cdot\mathbf{r}} \equiv \frac{1}{4\pi r_I} \]

Words. \(D_{00}\) = temporal gluon propagator of the lattice action (here Lüscher-Weisz), \(r_I\) = the improved distance, \(r_b\) = the bare distance \(|\mathbf{r}|\) in lattice units.

Derivation. In the continuum \(D_{00} = 1/\mathbf{k}^2\) and the Fourier transform is the Coulomb potential,

\[ \int\frac{d^3k}{(2\pi)^3}\frac{e^{i\mathbf{k}\cdot\mathbf{r}}}{\mathbf{k}^2} = \frac{1}{4\pi r}. \]

On the lattice \(\mathbf{k}^2\) is replaced by the lattice momentum \(\hat{\mathbf{k}}^2 = \frac{4}{a^2}\sum_i\sin^2\frac{ak_i}{2}\) (plus the improvement terms of the LW action) and the integration domain is the Brillouin zone. The resulting number is defined to be \(1/(4\pi r_I)\). Thus \(r_I\) is "the distance at which a continuum Coulomb potential would take the value the lattice actually produces". Replacing \(r\to r_I\) in the data is tree-level improvement; it removes the \(\mathcal{O}(g^0a^2)\) artefacts by construction and reduces the wiggles to the \(1\%\) level.

\[ K\left(\frac{r_b}{a},a,\beta_{\rm ref}\right) = \frac{V_S(r_I,a)-V_S(r_I,\beta_{\rm ref})}{V_S(r_I,\beta_{\rm ref})} \]

Words. \(\beta_{\rm ref}\) labels a finer reference lattice, which -- evaluated at \(r_b/a>\sqrt6\) where artefacts are negligible -- serves as a stand-in for the continuum answer. So \(K\) = relative deviation of the coarse lattice from the continuum, at fixed physical distance. Since \(V_S\) needs the additive renormalisation \(2C_Q\) and that is not known to sub-percent accuracy, the reference curve is shifted by a constant within its error until the ratio is \(\approx1\) in the safe region.

\[ K\left(\frac{r_b}{a}\right)=\sum_{i,j=1}^{i+j<5}b_{ij}\left(\frac{r_b}{a}\right)a^{2i}g_0^{2j} \]

Derivation of the allowed terms.

Only even powers of \(a\): Symanzik's analysis, and the fact that HISQ/LW has no \(\mathcal{O}(a)\) terms.
Powers of \(g_0^{2j}\) count loops.
Terms \(b_{i0}a^{2i}\) (i.e. \(j=0\), tree level) are omitted because the tree-level improvement \(r\to r_I\) already removed them.
Terms \(b_{0j}g_0^{2j}\) (i.e. \(i=0\)) are forbidden: they would survive \(a\to0\), but in the continuum rotational symmetry is exact and no \(r_b/a\)-dependent artefact can remain.

Hence both indices start at 1, and the truncation \(i+j<5\) is the practical order.

\[ V_{S,I} = \frac{V_S}{1+\langle K\rangle}, \qquad F_{S,I} = \frac{F_S}{1+\langle K\rangle} \]

Derivation. Eq. (B2) says \(V_S^{\text{lattice}} = V_S^{\text{cont}}\left(1+K\right)\). Inverting gives (B4); \(\langle K\rangle\) is the weighted average of \(K(\beta_{\rm ref})\) over several reference lattices, which reduces the statistical noise of any single choice. Eq. (B5) applies the same factor at \(T>0\): legitimate because the paper has shown that \(F_S \approx V_S\) for \(rT<0.3\), and the improvement is only ever used at \(r_b/a\lesssim4\).

\[ C^{\text{sub}}_{P,I} = \frac{1}{N_c^2}e^{-\frac{V_{S,I}-2F_Q}{T}}+\frac{N_c^2-1}{N_c^2}e^{\frac{V_{S,I}-2F_Q}{(N_c^2-1)T}} \]

This is exactly Eq. (30) with \(V_S\to V_{S,I}\) and with the Casimir-violating factors set to one (\(\delta_8=0\), \(\delta V_o=0\)), since we only want the change induced by the improvement.

\[ C_{P,I} = C_P+\frac{\langle K\rangle}{N_c^2-1}\frac{V_S}{T}\left[e^{-V_S/T}-C_P\right]+\mathcal{O}\left(\langle K\rangle^2\right) \]

Full derivation. Write \(A\equiv\frac{V_S-2F_Q}{T}\) and use \(V_{S,I}\simeq V_S(1-\langle K\rangle)\), so that \(\frac{V_{S,I}-2F_Q}{T} = A-\delta\) with \(\delta\equiv\langle K\rangle V_S/T\). Expand Eq. (B6) to first order in \(\delta\):

\[ C^{\text{sub}}_{P,I} = \frac{1}{N_c^2}e^{-A}(1+\delta)+\frac{N_c^2-1}{N_c^2}e^{\frac{A}{N_c^2-1}}\left(1-\frac{\delta}{N_c^2-1}\right) = C^{\text{sub}}_P+\delta\left[\frac{1}{N_c^2}e^{-A}-\frac{1}{N_c^2}e^{\frac{A}{N_c^2-1}}\right]. \]

Now eliminate the octet exponential using \(C^{\text{sub}}_P = \frac{1}{N_c^2}e^{-A}+\frac{N_c^2-1}{N_c^2}e^{\frac{A}{N_c^2-1}}\), i.e. \(\frac{1}{N_c^2}e^{\frac{A}{N_c^2-1}} = \frac{1}{N_c^2-1}\left(C^{\text{sub}}_P-\frac{1}{N_c^2}e^{-A}\right)\). Then

\[ \left[\ \cdot\ \right] = \frac{1}{N_c^2}e^{-A}-\frac{1}{N_c^2-1}\left(C^{\text{sub}}_P-\frac{1}{N_c^2}e^{-A}\right) =\frac{1}{N_c^2-1}\left[e^{-A}-C^{\text{sub}}_P\right], \]

where the last step used \(\frac{N_c^2-1}{N_c^2}+\frac{1}{N_c^2}=1\). Substituting back gives Eq. (B7). \(\square\)

\[ F^{\text{sub}}_{Q\bar Q,I} = -T\ln\left\{C^{\text{sub}}_P+\frac{\langle K\rangle}{N_c^2-1}\frac{F_S}{T}\left[C^{\text{sub}}_S-C^{\text{sub}}_P\right]\right\} \]

Derivation. Take (B7) and replace the two vacuum ingredients by their finite-\(T\) measured counterparts: \(V_S\to F_S\) and \(e^{-V_S/T}\to C^{\text{sub}}_S\) (which is \(e^{-F^{\text{sub}}_S/T}\) by Eq. 11). This is legitimate in exactly the range where the pNRQCD formula holds (\(rT<0.3\)), and it has the practical advantage that \(F_S\) and \(C_S^{\text{sub}}\) are available at every \(\beta\), whereas \(V_S\) is only available at the \(\beta\) values where \(T=0\) ensembles exist. Finally \(F^{\text{sub}} = -T\ln C^{\text{sub}}\). Valid to leading order in \(\langle K\rangle\).

The problem. The data live on an irregular grid: different \(N_\tau\) exist at different \(T\) and at different \(r/a\). To take a continuum limit at fixed \((rT,T)\) one first has to interpolate each \(N_\tau\) data set onto a common grid ("mock data"), then extrapolate in \(N_\tau\).

\[ F^{\text{sub}}(rT) = C+\frac{e^{-MrT}}{(rT)^\alpha}\sum_{n=n_{\min}}^{n_{\max}}A_n(rT)^n \]

Anatomy of the Ansatz.

\(\frac{1}{(rT)^\alpha}\): the leading short-distance power. \(\alpha=1\) for \(F^{\text{sub}}_S\) (Coulomb, \(1/r\)); \(\alpha=2\) for \(F^{\text{sub}}_{Q\bar Q}\) at short distance, dictated by the LO result Eq. (32) which goes like \(1/r^2\); \(\alpha=0\) for \(V_S-F_S\), which vanishes at \(r\to0\). So the exponents are not free -- they are read off from the weak-coupling results.
\(e^{-MrT}\): exponential damping, i.e. screening. \(M\) is the screening mass *in units of \(T\)*. It is forced positive, except that for \(T<235\) MeV fits with \(M=0\) are allowed too (there may be no resolvable screening in the confined phase).
\(\sum A_n(rT)^n\): a polynomial to absorb everything in between (running coupling, subleading terms). \(n_{\max}\le8\).
\(C\): a constant, allowing for the fact that \(F^{\text{sub}}\to0\) only in infinite volume; in a finite box the cancellation against \(2F_Q\) is incomplete.

The name says it: a Cornell potential \(a/r+\sigma r+C\), modified (arbitrary polynomial instead of \(a/r+\sigma r\)) and screened (the exponential).

\[ F^{\text{sub}}(rT) = \frac{e^{-M(T)rT}}{(rT)^\alpha}\sum_{m=m_{\min}}^{m_{\max}}\sum_{n=n_{\min}}^{n_{\max}}A_{n,m}(rT)^n\,T^m \]

Same shape, but the coefficients are now themselves polynomials (typically in \(T^{-3}\dots T^{0}\)) in the temperature, so that one fit describes the whole \((rT,T)\) surface at fixed \(N_\tau\). \(C\) is dropped, because no systematic \(T\) dependence was seen in the local \(C\)'s -- its effect is treated as a temperature-dependent uncertainty to be averaged out.

Why two schemes. The local fits are well constrained at short distance but ignore smoothness in \(T\); the global fits impose smoothness in \(T\) (hence much smaller errors at large \(r\)) but give little weight to the few short-distance points. Their difference is the paper's estimate of the systematic error.

\[ M(T) = M_0\left[1+M_1\ln\frac{T}{T_1}\right], \qquad T_1 = 1\ \text{GeV} \]

Reading. A screening mass in units of \(T\) should be almost constant, because \(m_D\propto g(T)T\) and \(g\) runs only logarithmically. So the natural model for \(M = m/T\) is a constant plus a logarithm; \(M_1\) is restricted to \((-1,+1)\) and comes out \(|M_1|\lesssim0.5\), an order of magnitude below \(M_0\) -- confirming that mild logarithmic running is enough. Normalising by \(T_1=1\) GeV keeps \(\ln(T/T_1)\in(-2,+2)\) over the studied range, so the coefficients stay natural.

\[ O(rT,T,N_\tau) = a(rT,T)+\frac{b(rT,T)}{N_\tau^2}+\frac{c(rT,T)}{N_\tau^4} \]

Derivation. The HISQ/tree action has leading cutoff effects \(\mathcal{O}(\alpha_s a^2,a^4)\) -- only even powers of \(a\). At fixed physical temperature, \(a = 1/(N_\tau T)\), hence \(a^2\propto1/N_\tau^2\) and \(a^4\propto1/N_\tau^4\). So \(a(rT,T)\) is the continuum value and \(b,c\) are the discretisation coefficients. The extrapolation is done independently at every point of the \((rT,T)\) plane, because cutoff effects vary strongly across regimes; the price is that the continuum surface is not guaranteed smooth.

\[ F^{\text{sub}}(rT,\dots) = C-\frac{AT}{rT}\exp\left[-MrT\right] \]

Derivation. At \(r\gg1/m_D\) the correlator is dominated by the exchange of the lightest single screened mode, so it must behave as a Yukawa: \(\propto e^{-mr}/r\). Writing \(m = MT\) and \(\frac{AT}{rT} = \frac{A}{r}\), this is exactly \(F^{\text{sub}} = C - A\frac{e^{-mr}}{r}\). The constant \(C\) is not physics: in infinite volume the correlator must approach \(L^2\) exactly, i.e. \(F^{\text{sub}}\to0\); in a finite box the cancellation is incomplete and leaves an offset. \(C\) is determined separately by fitting a constant to the noise-dominated tail (signal-to-noise \(<1\)), and subtracting it brings different volumes into agreement -- after which \(A\) and \(M\) are volume independent within errors.

Practice. Fit \(\ln\left[rT\left(F^{\text{sub}}-C\right)\right] = -\ln(AT)-MrT\) by linear regression, varying \((rT)_{\min}\) and \((rT)_{\max}\); \(M\) is read as the screening mass in units of \(T\), and a plateau in \(M\) as a function of \((rT)_{\min}\) is the signal that the asymptotic regime has been reached.

\[ C_{PL}(\beta,N_\tau,z) = \left\langle P(N_\tau,\mathbf{x})P^\dagger(N_\tau,\mathbf{x}+z\hat{e}_z)\right\rangle \]

with both Polyakov loops *integrated over the \(xy\) plane* before correlating. Since only the \(z\) separation remains, \(r = z\). Summing over a plane projects onto zero transverse momentum, which kills all excited states faster and improves the signal -- the same trick as zero-momentum projection in hadron spectroscopy.

\[ F^{PL}_{Q\bar Q}(rT,\dots) = C_0-A'T\exp\left[-MrT\right] \]

*Derivation (the \(1/r\) prefactor disappears).* Integrate the 3d Yukawa over the two transverse directions:

\[ \int d^2x_\perp\,\frac{e^{-m\sqrt{z^2+x_\perp^2}}}{\sqrt{z^2+x_\perp^2}} = 2\pi\int_0^\infty \rho\,d\rho\,\frac{e^{-m\sqrt{z^2+\rho^2}}}{\sqrt{z^2+\rho^2}} . \]

Substitute \(s = \sqrt{z^2+\rho^2}\), so \(s\,ds = \rho\,d\rho\) and \(\rho:0\to\infty\) becomes \(s:z\to\infty\):

\[ = 2\pi\int_z^\infty ds\,e^{-ms} = \frac{2\pi}{m}e^{-mz}. \]

So the plane-averaged correlator is a pure exponential with no power prefactor -- which is Eq. (D3), with \(M\) the same screening mass. The constant \(C_0\) now diverges in the continuum limit (the transverse integration is not normalised), but that is harmless since only the slope is wanted. Getting the same \(M\) from Eqs. (D1) and (D3) is a non-trivial consistency check that the asymptotic regime has really been reached.

Summary table of all equations

Eq.	What it is	Status
(1)	\(C_P = e^{-F_{Q\bar Q}/T}\)	definition / statistical mechanics
(2)(3)	lattice Polyakov loop from temporal links	discretisation
(4)	\(L_{\rm bare} = \langle P\rangle\)	definition
(5)	\(C_P^{\rm bare}=\langle PP^\dagger\rangle\)	definition
(6)	Coulomb-gauge singlet correlator	colour projection
(7)(8)(9)	divide by \(L^2\) to cancel the self-energy divergence	renormalisation
(10)(11)	\(F^{\rm sub}=-T\ln C^{\rm sub}\)	definition
(12)	\(L = e^{-F_Q/T}\)	statistical mechanics
(13)(14)	\(F = F^{\rm sub}+2F_Q\)	follows from (1),(8),(10),(12)
(15)	effective coupling from the force	definition, \(=\alpha_s\) at LO
(16)-(18)	short-distance thermal corrections, \(\mathcal{O}(g^4)\)	pNRQCD, quoted from Ref. [31]
(19)	soft/Debye correction, \(\mathcal{O}(g^5)\)	derivable from screened Coulomb, above
(20)	\(m_D\) at LO	one-loop \(\Pi_{00}(0,0)\)
(21)(22)	\(F^{\rm sub}_S\) at NLO in EQCD	quoted from Ref. [42]
(23)(24)	running-coupling correction, \(\delta Z_1\)	EQCD matching
(25)	\(m_D\) at NLO	quoted from Ref. [35]
(26)	singlet/octet decomposition of \(C_P\)	colour algebra, \(1\oplus8\)
(27)	\(L_A\)-weighted approximation of (26)	pNRQCD at \(r\ll1/T\)
(28)	octet piece alone	definition
(29)	\(L_A = L^{C_A/C_F}(1-\delta_8)^{C_A/C_F}\)	Casimir scaling
(30)	reconstruction of \(F^{\rm sub}_{Q\bar Q}\) from lattice inputs	derived from (27),(29),(13)
(31)	Casimir-scaling violation \(\delta V_o\)	\(\mathcal{O}(\alpha_s^3)\), Ref. [77]
(32)	\(F^{\rm sub}_{Q\bar Q}=-\frac{N_c^2-1}{8N_c^2}\frac{\alpha_s^2}{r^2T}\)	derived: \(1/r\) cancels, \(1/r^2\) survives
(33)	same, times \(e^{-2m_Dr}\)	two screened gluon exchanges
(A1)	pull for quark-mass dependence	statistics
(A2)	screening function \(S_1=-rF^{\rm sub}\)	linearises the Yukawa
(A3)	\(S_Q = -\partial F_Q/\partial T\)	thermodynamics
(B1)	improved distance \(r_I\)	lattice tree-level Fourier transform
(B2)(B3)	residual cutoff correction \(K\) and its interpolation	Symanzik counting
(B4)(B5)	\(V_{S,I}=V_S/(1+\langle K\rangle)\)	inversion of (B2)
(B6)(B7)(B8)	propagating the improvement to \(F_{Q\bar Q}\)	expansion of (30) in \(\langle K\rangle\)
(C1)(C2)(C3)	MSCP / GTDMSCP interpolating Ansätze	data-driven, shapes fixed by (32),(21)
(C4)	\(a+b/N_\tau^2+c/N_\tau^4\)	HISQ has only even powers of \(a\)
(D1)	asymptotic Yukawa fit	one-particle exchange
(D2)(D3)	plane-averaged correlator, pure exponential	transverse integration of the Yukawa

Everything measurable is a ratio of correlators, because that is what makes the linear divergence of the static Wilson lines cancel (Eqs. 7-9).
The Polyakov loop correlator is a colour average, not a potential. Its leading Coulomb term cancels between singlet and octet (Eq. 32), so \(F_{Q\bar Q}\) behaves as \(1/r^2\) and is exponentially harder to measure than \(F_S\).
Three regimes, three tools: \(r\ll1/T\) \(\to\) pNRQCD (Eqs. 16-19, 26-32); \(r\sim1/m_D\) \(\to\) EQCD (Eqs. 20-25, 33); \(r\gg1/m_D\) \(\to\) non-perturbative, only fits (Eqs. D1-D3).
Screening sets in at \(rT\approx0.3\) and is quantitatively perturbative (EQCD at NLO) for \(0.3\lesssim rT\lesssim0.6\), provided \(T\gtrsim300\) MeV.
Asymptotically, \(m_S\approx(1.6\text{-}2.0)\,m_D^{\rm NLO}\) and \(m_{Q\bar Q}\approx1.25\times2m_D^{\rm NLO}\): the pattern is perturbative, the normalisation is not.