Color screening in (2+1)-flavor QCD

We study correlation functions of spatially separated static quark-antiquark pairs in (2+1)-flavor QCD in order to investigate onset and nature of color screening at high temperatures. We perform lattice calculations in a wide temperature range, $140 \le T \le 5814\,{\rm MeV}$, using the highly improved staggered quark action and several lattice spacings to control discretization effects. By comparing at high temperatures our lattice results to weak-coupling calculations as well as to the zero temperature result for the energy of a static quark-antiquark pair, we observe that color screening sets in at $rT \approx 0.3$. Furthermore, we also observe that in the range $0.3 \lesssim r T \lesssim 0.6$ weak-coupling calculations in the framework of suitable effective field theories provide an adequate picture of color screening.


I. INTRODUCTION
As the temperature of strongly interacting matter increases to a pseudocritical temperature T c , a transition to a state with different properties than the vacuum at zero temperature occurs. Deconfinement of gluons and quarks, restoration of chiral symmetry, and screening of color charges are key properties of this thermal medium (cf. recent reviews, e.g., [1][2][3]).
For pure Yang-Mills theory [SUðN c Þ gauge theory], the Polyakov loop and the correlator of two Polyakov loops are the order parameters of the deconfinement transition, in which the ZðN c Þ center symmetry [4,5] of the pure Yang-Mills vacuum is broken. Since the center symmetry is broken by dynamical quarks already at zero temperature, neither the Polyakov loop nor the Polyakov loop correlator are order parameters of full QCD [6,7]. Nevertheless, the Polyakov loop correlator is a sensitive probe of deconfinement and color screening. In particular, it provides insight into the intricate screening properties of the thermal medium. The properly renormalized correlator of two Polyakov loops C P is related to the free energy of a static quark-antiquark (QQ) pair F QQ at separation r and temperature T [4,8], We choose to renormalize C P in such a way that F QQ for infinite separation matches the sum of the free energies of a decoupled static quark and antiquark using the renormalization of Ref. [9]. This is the same renormalization scheme that is usually used for matching the static energy at zero temperature for different lattice spacings.
At sufficiently high temperatures and short distances, the Polyakov loop correlators can be calculated using the weak-coupling expansion, i.e., the expansion in the gauge coupling g ¼ ffiffiffiffiffiffiffiffiffi ffi 4πα s p . In this regime, besides the scale 1=r and the Coulomb potential α s =r, there are three relevant thermodynamical scales that need to be considered: πT, gT (electric mass), and g 2 T (also called the magnetic mass). The leading order weak-coupling result for the Polyakov loop correlator has been known for a long time [4,28] and implies a non-Coulombic behavior of the static QQ free energy due to the cancellation of the singlet and octet contributions [4,28]. The next-to-leading-order calculation of the Polyakov loop correlator was performed by Nadkarni in the large distance regime, r ∼ 1=ðgTÞ [29], whereas the corresponding calculation in the short distance regime, rT ≪ 1, has been performed relatively recently [30]. In [31], one can find the state of the art expression for the Polyakov loop correlator at short distances, which is valid up to order g 7 and ðrTÞ 3 in the multipole expansion. In this regime, it turns out to be useful to organize the calculation by means of the effective field theory called potential nonrelativistic QCD (pNRQCD) [32][33][34]. For distances r ≳ 1=ðgTÞ, the behavior of the Polyakov loop correlator can be best understood in the framework of the dimensionally reduced effective field theory called electrostatic QCD (EQCD) [35] (see references therein for related works). At very large distances, r > 1=ðg 2 TÞ, the perturbative expansion breaks down due to the well-known Linde problem [36]. In Refs. [37,38], using EQCD, it was shown indeed that the very large distance behavior of the Polyakov loop correlator is affected by chromomagnetic screening and its correct asymptotic behavior is given by the lowest glueball mass at high temperatures [39,40]. The screening mass that governs the large distance behavior of the Polyakov loop correlator has been calculated in EQCD using lattice techniques and also investigated with weakcoupling methods in [41].
Most of the lattice QCD studies have focused on the large distance behavior of the Polyakov loop correlators and the extraction of the corresponding screening masses. The aim of this paper is to study the Polyakov loop correlator from small to large distances and see in what temperature and distance regime contact to weak-coupling calculations can be made. Such an analysis will also help to clarify at which distances color screening effects set in. As it was pointed out in the past, for understanding the temperature behavior of the Polyakov loop correlator, it is useful to consider the singlet correlator [8,30]. Therefore, in this paper, we will also present a detailed study of the singlet correlator and its comparison to weak-coupling calculations [31,42].
The question of the onset of color screening and applicability of weak-coupling expressions to the Polyakov loop and singlet correlators is not only of academic interest. There is an ongoing theoretical effort to understand quarkonium production in heavy-ion collisions (see, e.g., Refs. [43,44]), which requires the knowledge of in-medium quarkonium properties. Full QCD lattice determinations of in-medium quarkonium properties turn out to be difficult (see, e.g., Refs. [45,46]), which leads to the use of effective field theories [47][48][49][50][51][52][53]. Those, on the other hand, require one to know which scale hierarchy is relevant for the physically interesting case. The analysis presented in this paper aims at shedding some light on this issue.
The rest of the paper is organized as follows. In Sec. II, we cover technicalities of the lattice simulations. In Sec. III, we present our numerical results, the extraction of the continuum limit, and discuss some general features of the studied correlators, which do not require reference to explicit weak-coupling expressions. In Secs. IV and V, we conduct a systematic and quantitative comparison of the lattice results with the predictions from weak-coupling calculations. We finally use this knowledge for a discussion of the regime of asymptotic screening in Sec. VI. The last section of the paper contains our conclusions. Many technical details of the calculations are discussed in the Appendices. Some preliminary results from the present study have been reported in conference proceedings [54][55][56].

II. LATTICE QCD SETUP AND RENORMALIZATION
We perform calculations of the Polyakov loop correlator as well as of the singlet correlator of static QQ in (2 þ 1)flavor QCD at nonzero temperature on N 3 σ × N τ lattices with N τ ¼ 4, 6, 8, 10, 12, and 16 and aspect ratios of N σ =N τ ¼ 4 and 6 using the highly improved staggered quark (HISQ) action [57]. The calculations have been performed at the physical value of the strange quark mass m s and the light quark masses m l ¼ m s =20. The latter corresponds to a pion mass of m π ∼ 160 MeV in the continuum limit. We performed calculations in a wide beta range 5.9 ≤ β ¼ 10=g 2 0 ≤ 9.67, with g 0 being the lattice bare gauge coupling. The lattice spacing a has been fixed by the r 1 scale and we use the parametrization of r 1 =a given in Ref. [58]. Using this parametrization, we find that the above β range corresponds to a temperature range of 116 ≤ T ≤ 5814 MeV. The gauge configurations used in this study have been generated by the HotQCD Collaboration [7,58]. We also used the gauge configurations generated by the study of quark number susceptibilities at high temperatures [59,60]. Additional gauge configurations have been generated specifically for this study. Some of these gauge configurations have been used in the study of the renormalized Polyakov loop with the HISQ action [9,61]. A detailed account of the used gauge configurations, including the new gauge configurations by our (TUMQCD) collaboration is presented in Appendix A. Furthermore, in the high temperature region T > 350 MeV, additional gauge configurations have been generated with N τ ¼ 4, 6, 8, 10, 12, and 16, N σ =N τ ¼ 4 for β ¼ 7.03, 7.825, 8.0, 8.2, and 8.4, and light quark mass m l ¼ m s =5. This light quark mass corresponds to a pion mass of m π ≈ 320 MeV in the continuum limit. The aim of these calculations was to extend the continuum results on the correlators to higher temperatures. For temperatures T > 400 MeV, quark mass effects are expected to be negligible. We checked this explicitly in Appendix A. Finally, to study finite size effects for the large distance behavior of the correlators, we generated additional gauge configurations on 24 3 × 4 lattices (see Appendix A).
We have generated the new gauge ensembles using the rational hybrid Monte-Carlo (RHMC) algorithm and the MILC Collaboration code. Details on the HISQ action implementation in the MILC code can be found in Ref. [62].
In this paper, we study the correlation functions of the Polyakov loop. The (bare) Polyakov loop is defined in lattice QCD as the normalized trace of a temporal Wilson line WðN τ ; xÞ wrapping around the time direction once, where U 0 ðx ¼ ðτ; xÞÞ are the temporal link variables. We denote the expectation value of the bare Polyakov loop averaged over the spatial lattice volume as The bare Polyakov loop correlator is defined as the ensemble average (spatial average over x included) C bare P ðβ; N τ ; rÞ ¼ hPðN τ ; xÞP † ðN τ ; x þ rÞi: The superscript "bare" refers to the fact that the quantity is not renormalized. In addition, we study the singlet correlator in the Coulomb gauge [8,14], defined as C bare S ðβ; N τ ; rÞ ¼ as well as the cyclic Wilson loops [15,63,64] W bare S ðβ; N τ ; rÞ where Sðτ=a; x; rÞ is a spatial Wilson line between the points ðτ; xÞ and ðτ; x þ rÞ. These quantities can be considered as meson correlation functions composed of a static quark and antiquark separated by a distance r [15,63] and evaluated at Euclidean time τ ¼ 1=T [15,63]. The subscript bare refers again to the fact that these are unrenormalized quantities. At very large distances, C bare P ðβ; N τ ; rÞ, C bare S ðβ; N τ ; rÞ, and W bare S ðβ; N τ ; rÞ approach ½L bare ðβ; N τ Þ 2 in the infinite volume limit [8,15,65]. 1 These bare quantities contain ultraviolet (UV) divergences that are regularized on the lattice. The continuum limit requires looking for combinations of bare quantities that are free of UV divergences.
The renormalization of C bare P ðβ; N τ ; rÞ is simple: the only divergences here are the ones associated with the two Polyakov loops. Therefore, the normalized Polyakov loop correlator C sub P ðβ; N τ ; rÞ ¼ is free of UV divergences and the continuum limit can be taken. In particular, this implies that Eq. (7) can be written directly in terms of renormalized quantities, namely, C sub P ðβ; N τ ; rÞ ¼ where C P is the renormalized Polyakov loop correlator and L is the renormalized Polyakov loop. The normalized Polyakov loop correlator C sub P is simply the connected part of the Polyakov loop correlator and contains the complete information about the in-medium modification of the QQ interaction. The disconnected part of the renormalized correlator, i.e., L 2 , contains information about infinitely separated static quarks or antiquarks that interact exclusively with the medium, but not with each other.
A similar observation holds for the singlet correlator and we can use the normalized correlation function also in this case, C sub S ðβ; N τ ; rÞ ¼ C bare S ðβ; N τ ; rÞ ½L bare ðβ; N τ Þ 2 ¼ C S ðβ; N τ ; rÞ ½Lðβ; N τ Þ 2 : ð9Þ The above expression, which is defined in the Coulomb gauge, is free of divergences and has a well-defined continuum limit. Such a statement does not hold for general gauges: in covariant gauges, the normalized correlator will be UV divergent [42]. Another reason for using the Coulomb gauge is that weak-coupling calculations are available in this case [30,31,42].
Since the normalized Polyakov loop and singlet correlators are ultraviolet finite, we can define the subtracted free energy of a static QQ pair and the so-called (subtracted) singlet free energy by taking the logarithm of these correlators, 1 Here we take into account that in QCD L is real. F sub QQ ðr; T; aÞ ¼ −T ln C sub P ðβ; N τ ; rÞ; ð10Þ F sub S ðr; T; aÞ ¼ −T ln C sub S ðβ; N τ ; rÞ: ð11Þ Here we have traded β and N τ for the physical temperature T ¼ 1=ðN τ aÞ and the lattice spacing a. Up to a finite constant, Eqs. (1) and (10) define the same free energy. From Eq. (8), it follows that this constant is given by −2T ln L. The expectation value of the renormalized Polyakov loop L is related to the free energy of an isolated static quark F Q at temperature T, The constant −2T ln L is therefore the free energy of two static quarks or antiquarks that are at infinite separation and, thus, can only interact with the medium. The regularized UV divergence of L bare is exactly the square root of the regularized UV divergence of C bare P or C bare S . Hence, by adding twice the free energy of an isolated static quark F Q , determined in Ref. [9], we get the free energy of a static QQ pair as well as the singlet free energy, Equation (13) trivially follows from Eqs. (1), (8), (10), and (12). Equation (14) is the analog for the singlet free energy. Thus, the problem of renormalizing the Polyakov loop and singlet correlators is reduced to the problem of renormalizing L bare . For the latter, we use the procedure described in Ref. [9] and update it with the most recent lattice results. The updated renormalization of L bare is presented in Appendix A. The above procedure also turns out to be very convenient when performing the continuum extrapolations discussed in the next section.
Another way to study in-medium modifications of the interaction of static quark and antiquark is to consider cyclic Wilson loops [15,63,65,66]. The renormalization of the cyclic Wilson loops, however, is more complicated. In this case, in addition to the self-energy divergences present in the Polyakov loop correlators, there are also self-energy divergences associated with the spatial Wilson lines SðN τ ; x; rÞ and S † ð0; x; rÞ, as well as intersection divergences [64]. On the lattice, there are also cusp divergences for off-axis separations. Rectangular Wilson loops evaluated for τ < 1=T have cusp divergences at the four corners of the loop. These cusp divergences are absent for cyclic Wilson loops (i.e., for τ ¼ 1=T) in the continuum, because the corresponding contour is smooth. However, on the lattice, one can also consider cyclic Wilson loops for noninteger multiples of the lattice spacing by using Wilson lines that are not aligned with any of the spatial lattice axes. Such off-axis spatial separations, e.g., r=a ¼ ffiffi ffi 2 p , involve cusps. In other words, for off-axis separation, the cyclic Wilson loops have additional cusp divergences, while for on-axis separation, such divergences are absent. To simplify the analysis, we only consider onaxis separation. Constructing the spatial part of the cyclic Wilson loops from the smeared gauge links can reduce the size of the UV divergences. Therefore, we calculate the cyclic Wilson loops with spatial lines obtained from smeared gauge links. We apply several steps of hypercubic (HYP) smearing [67] and use the same parameters for HYP smearing as in Ref. [67].

III. LATTICE RESULTS
In this section, we discuss the dependence of the Polyakov loop correlator, defined in Eq. (5) and the singlet correlator, given in Eq. (6), on the separation of the static quark-antiquark pair at different temperatures. The lattice results for the Polyakov loop correlator and the singlet correlator are functions of three variables: the distance r, the lattice spacing a, and the temporal extent N τ . The temperature T is related to N τ and a through aN τ ¼ 1=T; by trading a for the temperature T, we parametrize cutoff effects by the N τ dependence of the correlators at fixed temperature. In total, we analyze a full set of more than 10 000 data (expectation values) from more than 150 gauge ensembles. The temperature range covers the confined phase, the crossover region, and a strongly as well as weakly coupled quark-gluon plasma phase. The separation range covers asymptotic freedom at short and color screening at large distances.

A. Outline of the analysis
The correlator data for different N τ are available only for heterogeneous sets of temperatures and for heterogeneous sets of separations at each temperature. Therefore, the comparison of correlators obtained for different N τ and eventually the continuum extrapolation requires either a three-dimensional, fully global fit or, in one way or another, an interpolation in the separation r (or in rT) and in the temperature T using smoothing functions in order to obtain an intermediate, homogeneously spaced grid of mock data. We perform separate fits in the separation and the temperature sequentially (local fits) or simultaneous fits in the separation and the temperature (global fits). Smoothness in both separation and temperature is a physical constraint in the absence of sharp phase transitions, but inevitably leads to a shrinkage of the errors in the process. Hence, we have to expect that errors of the mock data may be underestimated.
Since no analytic results are available that cover the whole range spanned by the available lattice data, any interpolation introduces some degree of bias. In order to minimize this bias, we identify classes of apparently sensible model functions. For small separations, such model functions are constrained by asymptotic freedom and expectations based on the leading-order weak-coupling results. For larger separations, such model functions are constrained by the dominance of massive one-particle exchange processes, i.e., exponential damping due to color screening. We do not enforce exponential damping for lower temperatures, arguing that confinement may still be unabated or that we may fail to resolve the damping due to a lack of data with sufficient accuracy at large enough separations. For even larger separations, subtracted free energies must approach zero due to cluster decomposition, unless statistical fluctuations (for smaller ensembles) or artifacts of the finite volume lead to an incomplete cancellation and, thus, to a constant value. We account for all of these effects in the interpolations discussed in Appendix C and use a data-driven procedure to identify the most adequate (i.e., most simple) interpolating functions. We note that in the case of F sub QQ we actually find that the lattice data successfully override our inherent bias at low temperatures. Namely, we identify Coulombic interactions despite our explicit bias towards non-Coulombic interactions (see the discussion in Sec. III C).
As ensembles at different temperature T (or at different N τ ) are statistically independent, correlated fluctuations that affect the size of error estimates are relevant only for an interpolation in the separation. Since we are interested in the behavior of the correlators in a wide range of distances r, down to the shortest ones, it is important to take care of the lattice artifacts present at short distances before performing the interpolation in r. To do this, we utilize the results of the lattice calculation of the static QQ energy at zero temperature V S (published in Refs. [58,68] and the procedure described in Ref. [69]). We give a detailed account of the underlying procedure in Appendix B. When discussing the T ¼ 0 static energy throughout this section, we always use these nonperturbatively improved lattice results. This improvement procedure, however, introduces systematical uncertainties at short distances that we have to estimate heuristically. Such systematical uncertainties are manifestly uncorrelated and orders of magnitude larger than the corresponding statistical errors. That is why the total uncertainties of the data underlying the interpolations (including short distances) are always dominated by uncorrelated errors. Hence, uncorrelated fits aiming for χ 2 =d:o:f: ≃ 1 are an appropriate procedure for analyzing these data, where d.o.f. represents degrees of freedom.
We generalize the improvement procedure to correlators at T > 0. Here, in particular, for ensembles with smaller N τ , we have to be aware that medium effects may become important in the short distance range, where we apply the nonperturbative improvement. For this reason, we have to expect that we may very well underestimate the associated uncertainties and account for this by enlarging the errors if we cannot interpolate smoothly between data at short distances with χ 2 =d:o:f: ≃ 1. After correcting for short distance lattice artifacts, we perform interpolations of rF sub S and r 2 TF sub QQ in rT and T using splines or polynomial fits with exponentially decaying prefactor of the form exp ½−MrT, where M is a fit parameter. We perform uncorrelated fits on the individual jackknife bins and determine the errors of the fit from the respective variance of the fit model. These interpolations are only weakly constrained for the largest separations.
In a secondary interpolation at fixed rT, we use, in the absence of justified model assumptions, smoothing splines to interpolate the mock data in T using uncorrelated fits, since ensembles at different T are statistically independent. Correlators at lower temperatures usually have consistently larger uncertainties than correlators at higher temperatures. For this reason, we have to separately fit in overlapping low and high temperature ranges in order to be sensitive to the lower temperatures. As we typically have about 30 different temperatures for each N τ , we find that smoothness under variation of the temperature can provide tighter constraints for the rT dependence than the actual correlations in each ensemble at a fixed temperature. This justifies the previous use of uncorrelated fits. In order to directly incorporate the smoothness under variation of T in the analysis of the rT dependence, we also perform global fits. Namely, we fit the rT and T dependence simultaneously in separate but overlapping low and high temperature intervals for each N τ . In these global fits, we model the temperature dependence heuristically with (inverse) monomials in T and with a T-dependent Ansatz for the screening mass. The global fits yield much smaller uncertainties at larger distances than the local fits, but are less effectively constrained by the data at short distances. Namely, the few data at short distances with large systematic errors have only minor weights in the global fits. Eventually, we estimate systematic uncertainties from the differences between local and global fits, since their strengths and deficiencies are complementary.
Using the mock data for these quantities at the intermediate values of rT and T obtained from interpolations, we perform fits of the form a þ b=N 2 τ þ c=N 4 τ to obtain the continuum limit. This form is appropriate for the HISQ action, which has cutoff effects only at even powers of the lattice spacing. We encounter cases where the expected scaling behavior is observed on a set of coarser lattices, i.e., smaller N τ , but apparently not realized for finer lattices. We understand that this may be caused by underestimating the errors of the mock data and rescale the errors in order to obtain a more realistic (i.e., larger) uncertainty of the continuum result in these cases. As the lattice data cover vastly different physical regimes in the ðrT; TÞ plane, it is unsurprising that the cutoff effects vary significantly over the whole data set. In the absence of prior knowledge 2 about the variation of cutoff effects, we conservatively perform a separate fit for each point in the ðrT; TÞ plane. The coefficients b and c in these fits sometimes were set to zero and the range in N τ was restricted to estimate the systematic errors of the extrapolations. These systematic errors can be accidentally large for individual points. However, only subsets of the N τ range permit access to the various regions in the ðrT; TÞ plane. As a consequence, we have to combine different continuum extrapolations in various regions and, thus, obtain minor discontinuities in the continuum result. The details of this analysis are discussed in Appendix C. The procedure works as long as the statistical errors on F sub S and F sub QQ are sufficiently small. At large distances, F sub S and F sub QQ approach zero and the relative statistical errors become very large. For these distances, our procedure does not work and we cannot provide a continuum result for the free energy and the singlet free energy. This limits the range of rT where continuum results can be provided. Finally, to obtain the free energy of a static QQ pair and the singlet free energy, we add the continuum limit of 2F Q to the continuum results of F sub S and F sub QQ . We note that for most of the range the overall error of F S and F QQ is dominated by the uncertainty of their asymptotic value 2F Q .

B. Continuum results
Our continuum results for F S are shown in Fig. 1 (left) and compared to the zero temperature static QQ energy V S , shown as black pentagons connected by a thin black band. These represent the nonperturbatively improved static energy for five different lattice spacings in overlapping intervals. The details of the composition of these data are discussed in Appendix C. The numerical data for the zero temperature static energy are taken from Refs. [7,58,68]. The asymptotic values, i.e., 2F Q , are shown as horizontal bands. As one can see from the figure, at high temperatures, T > 200 MeV, continuum results are available up to sufficiently large distances to see the approach to 2F Q . At short distances, the singlet free energy is temperature independent and agrees with the vacuum static QQ energy. At larger distances, we see the onset of medium effects. As the temperature is increased, the onset of the medium effects happens at shorter and shorter distances, but there is still an extended region in r where F S agrees with the zero temperature static energy. This is expected from weakcoupling calculations in pNRQCD [31].
The continuum results for F QQ are shown in Fig. 1  (right). Again, only for T > 200 MeV our continuum extrapolation extends to large enough distances to make connection to the asymptotic value of 2F Q . To compare to the zero temperature static energy, we subtract the constant T ln N 2 c ¼ T ln 9, which is due to the standard normalization convention of the Polyakov loop correlator defined by Eqs. (2) and (5), since usually the Polyakov loop correlator is defined without the 1=N 2 c ¼ 1=9 normalization of the trace (see, e.g., [63]). For the T ¼ 0 static energy, we use the same composite band as in the figure for F S . We extend this band to even shorter distances using the continuum limit of F S at very short distances as an estimate of V S . We show this estimate with a light gray color to keep the distinction clear and discuss the details of this band in Appendix C. For F QQ , we see a much larger temperature dependence already at relatively short distances and for temperatures T > 800 MeV no clear agreement with the vacuum static QQ energy can be seen. Qualitatively this behavior of the free energy and singlet free energy of QQ is very similar to the results of previous calculations performed at fixed lattice spacing [21][22][23][24][25][26]. The reason for this behavior of the QQ free energy is related to the presence of a color octet contribution and will be discussed in Sec. V. Finally, we mention that our continuum results for F QQ agree well with the continuum extrapolated results presented in Ref. [26] for T ≤ 400 MeV. Our analysis, however, extends to higher temperatures.

C. The effective coupling
At T ¼ 0 it is customary to study the force between the static quark and antiquark defined in terms of the static QQ energy V S ðrÞ or, equivalently, the corresponding effective coupling where C F ¼ 4=3 is the Casimir operator in the fundamental representation of SU(3). Following Ref. [70] we can generalize this approach to the free energy and singlet free energy and obtain the corresponding effective couplings. The reasoning behind this generalization is that, as we see in Fig. 1, F S and F QQ − T ln 9 are similar to V S at sufficiently small distances. At such distances, a sensible definition of an effective coupling seems possible. First, we discuss our numerical results for α QQ ðrÞ from the singlet free energy, which are shown in Fig. 2. The effective coupling defined in terms of F S is temperature independent at short distances and agrees with the result obtained from a numerical derivative of the static energy V S (details are given in Appendix C). At short distances, the running of the coupling is controlled by the scale 1=r. The effective coupling reaches a maximum at some distance r ¼ r max ðTÞ and then decreases, indicating the onset of color screening. Furthermore, it turns out that r max approximately scales like 0.4=T. The distance r max roughly separates the region that is dominated by vacuum physics from the one where screening sets in. 3 This will be addressed more quantitatively in the next section. It is interesting to see that the value of α QQ ðr max ; TÞ is quite large for T < 320 MeV, implying that the physics is strongly coupled. At higher temperatures, α QQ ðr max ; TÞ < 0.5, indicated by a horizontal line in the figure, so that weakcoupling methods might be applicable.
In the case of the free energy F QQ , we expect at short distances and for weak coupling two possible regimes. When T ≪ α s =r, F QQ should exhibit a typical Coulombic behavior 1=r, whereas for T ≫ α s =r, due to a cancellation between color singlet and color octet contributions, the behavior of F QQ as a function of the distance should be like 1=r 2 . Indeed, the subtracted free energy, defined in Eq. (10), multiplied by −r 2 T, shown in Fig. 3, 4 seems to exhibit both regimes. At short distances and for T < 500 MeV, −r 2 TF sub QQ raises with the distance r, a behavior consistent with F sub QQ ∼ 1=r, while for T > 500 MeV, the available data for −r 2 TF sub QQ appear less sensitive to r, a behavior consistent with F sub QQ ∼ 1=r 2 . Analytical results and a detailed quantitative analysis can be found in Sec. VA. Here our purpose is to justify the use of Eq. (15), at short distances and for T ≲ 500 MeV, to define an effective coupling also in the case of the free energy F QQ . Eventually, also this effective coupling can be compared with the coupling obtained from the QQ static energy at zero temperature.
The statistical errors on our lattice results for F QQ are much larger than in the F S case and therefore we calculate the corresponding effective coupling only for finite N τ . In Fig. 4, we show results for N τ ¼ 12 and 16. Following the above discussion, we look at temperatures T ≲ 500 MeV and at distances short enough that −r 2 TF sub QQ raises with the distance. Because of lack of accuracy of the derivative at the first data point for each temperature, we have to exclude where finite temperature effects are smaller than the statistical errors. The horizontal bands correspond to the r → ∞ limit of F S and F QQ , i.e., 2F Q . The subtraction of T ln 9 in the right panel is due to the normalization convention used for F QQ as discussed in the text. derivatives for rT smaller than the center of the interval between the first two data points for each temperature. The points in Fig. 4 satisfy the above requirements. The figure supports a description where, in the considered temperature range and for the shortest distances, the free energy of a QQ pair is dominated by the singlet contribution and agrees well with the T ¼ 0 static energy.

D. Cyclic Wilson loops
In this subsection, we present our numerical results for cyclic Wilson loops W bare S . Cyclic Wilson loops being gauge invariant offer an attractive possibility to study the medium modifications of QQ interactions. From previous lattice results, we know that W bare S approaches ½L bare 2 at large distances [6,65]. This is also expected based on weakcoupling calculations [42]. Like the Coulomb gauge singlet correlator, W bare S is also dominated by the singlet contribution at small distances [64,66]. Moreover, weak-coupling calculations show that W bare S is screened at large distances like the Coulomb gauge singlet correlator. However, the study of cyclic Wilson loops is complicated by the presence of UV divergences. In the continuum theory, the renormalization of cyclic Wilson loops has been worked out for short distances, rT ≪ 1 [64,66]. However, on the lattice, we do not have enough data points at truly short distances to test these ideas. In part, this is due to the fact that off-axis separation cannot be used in the analysis (see the discussion in Sec. II). Therefore, we attempt to compare the unrenormalized smeared cyclic Wilson loops with the Coulomb gauge correlator. The smearing of the spatial links in W bare S is expected to reduce the size of UV divergences due to the self-energy, which are the largest divergences that affect the spatial lines. As we will see, for smeared cyclic Wilson loops, the UV divergences are sufficiently small at the considered values of N τ to allow for a meaningful comparison.
We discuss our numerical results for cyclic Wilson loops in terms of F sub W ¼ −T lnðW bare S =½L bare 2 Þ. The analogy with Eq. (11) is obvious. We show the logarithm of −rF sub W in Fig. 5. There, we also show the usual Coulomb gauge singlet correlator in direct comparison. It is evident that −rF sub W from the unsmeared Wilson loop decays much more rapidly than the singlet correlator in the Coulomb gauge. This is due to the self-energy divergence of the spatial lines of the loop, which effectively increases the screening mass. As we increase the number of HYP smearing iterations, the size of the self-energy divergence in the spatial line is gradually reduced. After a few steps of HYP smearing, we obtain results for −rF sub W that agree with the results from the singlet correlator in the Coulomb gauge. The number of HYP smearing that is needed to sufficiently reduce the size of the self-energy divergence is larger for smaller lattice spacing. For higher temperatures, and thus finer lattice spacings at fixed N τ , more smearing iterations are required in general, which can be seen by comparing the lower temperature (T ≈ 200 MeV, upper panel in Fig. 5) with the higher temperature (T ≈ 480 MeV, lower panel). We see that larger N τ (finer lattices) requires larger numbers of smearing iterations to achieve a comparable effect. These findings are in qualitative agreement with the study performed in the SU(2) gauge theory [15].
Though both correlators are affected by different technical difficulties, the accessible information regarding the QQ interaction and the properties of the quark-gluon plasma seem to be similar.

IV. COMPARISON OF THE SINGLET FREE ENERGY TO WEAK COUPLING
In this section, we compare the lattice results for the singlet free energy with the predictions of the weakcoupling calculations. These comparisons are the key for a deeper understanding of the mechanisms behind the onset of thermal effects with increasing separation of the QQ pair and provide a tool to distinguish between electric and magnetic screening on the lattice.
The numerical results discussed in the previous section suggest that medium effects for the singlet free energy are small for rT < 0.4. Furthermore, weak-coupling calculations have been worked out up to order g 5 for rT ≪ 1 [31] and to next-to-leading order (NLO) for r ∼ 1=m D [42], which is accurate up to order g 4 [31]. Hence, we will compare to the weak-coupling results separately for short (rT ≪ 1) and large distances (rm D ∼ 1). This separation of different regimes is, strictly speaking, valid only for asymptotically small couplings g ≪ 1. However, as we will see below, the scale separation and the corresponding effective field theory calculations can also help to understand the lattice results when g ∼ 1.
A. Thermal corrections for small r Here we compare the lattice results with the weakcoupling calculations for rT ≪ 1. In this case, the singlet free energy is given by the energy of a static QQ pair in the vacuum V S , plus thermal corrections, which, however, are expected to be small. For this reason, it makes sense to perform the comparison with the lattice results in terms of the difference V S − F S . The comparison of the lattice and weak-coupling results in terms of this difference has several advantages. First, one does not have to deal with the complications associated with the perturbative calculation of V S , among others, infrared logarithms (see Refs. [69,71]). Second, lattice artifacts largely cancel out in this difference, making the continuum extrapolations of this quantity simpler. Finally, the ultraviolet additive renormalization does not enter in V S − F S . We show our continuum results for V S − F S in Fig. 6. The details of the continuum extrapolations are discussed in Appendix C. The difference V S − F S in physical units is small at short distances and roughly temperature independent. As the temperature increases, the distances where medium effects are significant become smaller. Since we  expect that thermal corrections should be proportional to the temperature, we present our results in temperature units in the right panel of Fig. 6. We see that for T > 300 MeV the difference V S − F S in temperature units approaches at short distances a constant. As we will explain in the next paragraph, this behavior is expected based on the weakcoupling calculations. Now we discuss the comparison of V S − F S with the weak-coupling calculations. From the point of view of weak-coupling calculations, short distances imply the scale hierarchy 1=r ≫ T. In this regime, pNRQCD can be used. Moreover, the above scale hierarchy allows two further possible scale hierarchies: α s =r ≫ T and T ≫ α s =r. We discuss the second case, since here explicit calculations for V S − F S are available up to unknown corrections of order g 4 ðrTÞ 5 and g 6 [31]. The result reads with contributions from the nonstatic gluons ΔF g , from the quarks (N f massless quark flavors) ΔF f , and from the soft (Debye mass) scale ΔF S , ; ð17Þ It is clear from the above expression that V S − F S vanishes in the r → 0 limit. This seems to be in apparent contradiction with the lattice results shown in Fig. 6. However, as remarked in Ref. [31] (see also the singlet normalization in [33]), there is a possible order g 6 contribution, which is proportional to the temperature and independent of r. Such a contribution would naturally explain the trend seen in the lattice results in Fig. 6 (right). Note that the ðrTÞ 2 term in Eq. (19) has the opposite sign with respect to the ðrTÞ 2 terms in Eqs. (17) and (18), leading to a partial cancellation. This partial cancellation is in fact responsible for relatively small medium effects in the considered distance region. In Fig. 7, we compare the lattice QCD results for T ¼ 500 MeV to the perturbative calculation reproduced in Eq. (16). We use Λ MS ¼ 320 MeV and the two-loop running coupling for this comparison. This choice of ΛM S and running coupling gives α  the lattice QCD determination of α N f ¼3 s based on the static energy [69], as well as to the lattice determination based on the charmonium correlators [72]. For the renormalization scale, we choose μ ¼ 2πT and add a small constant to the above weak-coupling results, which mimics the possibly missing g 6 contribution such that the weak-coupling calculation agrees with the lattice at the shortest distance. We vary the renormalization scale by a factor of 2 around the central value to estimate the uncertainty of the perturbative result. This uncertainty is shown as a band in Fig. 7. Because of the partial cancellations between ΔF g þ ΔF f and ΔF S , the net temperature shift of V S − F S is relatively small for rT ≲ 0.3 in agreement with the lattice results. Thus weak-coupling results obtained for rT ≪ 1 enable us to explain the lattice findings on F S for small r, i.e., for r ≲ 0.3=T, in particular, their small temperature dependence. For larger distances, ΔF S becomes numerically the dominant correction, even though it is parametrically small (order g 5 ) compared to ΔF g and ΔF f , which are of order g 4 . This means that screening effects become important, and the hierarchy of scales 1=r ≫ T ≫ m D is not fulfilled anymore. We need to consider a weak-coupling expansion that does not assume such a hierarchy and resums screening effects right from the start. This is discussed in the next subsection.
There is no perturbative calculation of F S − V S for the case α s =r ≫ T. However, from the power counting, we expect thermal corrections to be small [48]. Furthermore, we also expect an r-independent term proportional to T appearing at order g 6 , since such a term possibly arises from the matching of pNRQCD to nonrelativistic QCD [31]. Therefore, several qualitative features of F S − V S discussed above should also hold in this hierarchy of scales. This is quite different from the case of F QQ , where the two regimes α s =r ≪ T and α s =r ≫ T differ at leading order in the functional dependence on r.

B. Electric screening of the singlet free energy
In this subsection, we compare the lattice results for the singlet free energy with the weak-coupling prediction in the regime r ∼ 1=m D , i.e., in the electric screening regime. As mentioned in the Introduction, in this case, EQCD can be used and the correlator can be calculated perturbatively. The physics in this regime is controlled to a large extent by the mass parameter of EQCD, m D , which is often called the Debye mass. At leading order, The comparison of the lattice and the weak-coupling results is performed in terms of the subtracted free energy F sub S defined in Eq. (11). The comparison in terms of F sub S is most convenient, since F S is very close to 2F Q in this distance regime, and the weak-coupling calculations of 2F Q are difficult [73]. Our continuum results for F sub S are shown in Fig. 8. We clearly see that at rT > 0.25 F sub S decreases with increasing r. This fact is consistent with the observation made in the previous section that medium effects become large at rT > 0.3. The subtracted singlet free energy in the Coulomb gauge has been calculated at NLO in Ref. [42] using EQCD. We use Eq. (3.22) from that reference and omit the two-gluon exchange contribution proportional to g 4 =r 2 , which is of order g 6 in the regime r ∼ 1=m D [31]. The result reads where E 1 ðzÞ is the exponential integral, i.e., There are two corrections to the above result that are formally of higher order but numerically important. The first one is the correction due to the running coupling [31,42], which can be written in the form [31] δF sub S ðr; TÞ ¼ − with The other correction is the correction to m D . The NLO result for the Debye mass, which has been calculated in Ref. [35], reads In Fig. 8, we compare the weak-coupling results at LO and NLO with the continuum extrapolated lattice results. For the LO one, we use F sub S j LO þ δF sub S , and for the NLO one, we use F sub S j NLO þ δF sub S . For m D , we always use the above NLO result. As before, we use Λ MS ¼ 320 MeV and the two-loop running coupling. We vary the scale μ from πT to 4πT to account for the uncertainty of the weakcoupling result. First of all, we see that the scale dependence is much reduced for higher temperatures. For T ≳ 600 MeV, it is sufficiently small to allow for an accurate comparison with the lattice. Moreover, the central value of the NLO result is always higher than the central value of the LO result. Within the uncertainty band, due to the variation of the scale, we see good agreement between the NLO result and the lattice data up to at least rT ≈ 0.6. For rT ∼ 0.3-0.6, the shape of the lattice results and the weak-coupling results is quite similar. At larger distances, however, the falloff of the lattice results is faster than the one of the NLO prediction. This is most likely due to the fact that, for these distances, the physics is affected by static chromomagnetic fields and, thus, it is nonperturbative. Based on previous calculations performed in the pure gauge theory, we expect that the effective screening mass governing the large distance behavior of the singlet correlator becomes significantly larger than m D due to chromomagnetic effects [14,70,74,75]. This would explain the change in the slope at large distances that we observe. The behavior of F S in the large distance regime will be discussed in Sec. VI. However, it is important to stress that chromoelectric screening effects are important already at distances rT < 0.6.

V. COMPARISON OF THE POLYAKOV LOOP CORRELATOR TO WEAK COUPLING
In this section, we compare the lattice results for the Polyakov loop correlator C P or, equivalently, the free energy of a static QQ pair F QQ , with the predictions of the weak-coupling calculations. Here and in the following, we use small letters for the indices to mark color singlet or color octet objects defined in the framework of weakcoupling calculations and capital letters for the indices to mark color singlet or color octet objects defined in the context of nonperturbative lattice calculations. As in the case of the singlet free energy, we will perform this comparison separately for the short distance regime and for the screening regime, where r ∼ 1=m D .

A. The Polyakov loop correlator at short distances
We start by assuming the case 1=r ≫ T ≫ m D . As it was discussed in Sec. III, the free energy of a static QQ pair, F QQ − T ln 9, shows large deviations from the zero temperature static energy already at short distances. One way to understand this is in the language of pNRQCD. In the pNRQCD framework, the Polyakov loop correlator at short distances can be written up to corrections of order α 3 s ðrTÞ 4 as [31] where f s and f o are gauge-invariant singlet and octet free energies defined in terms of pNRQCD correlators [see Eqs. (71) and (72) of Ref. [31]]. A good approximation of (26) is respectively, the quark-antiquark color singlet and octet potentials at T ¼ 0, and L A is the adjoint Polyakov loop expectation value. Equation (27) is valid up to corrections of order α 3 s [cf. Eqs. (79) and (80) in Ref. [30]]. From these equations, it is clear that the temperature dependence of the Polyakov loop correlator, and thus of the QQ free energy, comes, through the color octet contribution, mostly from L A .
To test to what extent the above formulas work in the temperature and distance ranges considered in this study, we define the octet contribution to the Polyakov loop correlator as C o ≡ C P − expð−V s =TÞ=N 2 c . In Fig. 9 (left panel), we show C O (normalized by 8C P =9), which is a proxy for C o directly defined on the lattice. Specifically, this means that we are using lattice data for C P (on N τ ¼ 8 lattices) and the zero temperature static energy V S , calculated on the lattice in Refs. [58,68], as a proxy for V s . In this way, we avoid the complications associated with the perturbative calculation of V s . We see that C O increases with increasing temperatures and increasing distances. This can be understood by considering that V o is positive. Furthermore, at low temperatures, L A is small, thus leading to additional suppression of the octet contribution, while at high temperatures, L A is of order one [76]. According to Eq. (26), we can also parametrize C o as Thus, the proxy C O defines a proxy F O for the octet free energy f o . The results for F O are also shown in Fig. 9 (right panel). We see that F O decreases with increasing temperatures. Moreover, for distances r < 0.1 fm, we see hints for the presence of the repulsive tail V o .
To validate the pNRQCD framework on a quantitative level, we reconstruct the subtracted free energy of a static QQ pair using the pNRQCD decomposition formula (27) and lattice results for the zero temperature static QQ energy V S , taken from Refs. [58,68], as a proxy for V s , and for the Polyakov loop expectation value in the fundamental representation, L F ≡ L, taken from Ref. [9]. First, we write for the adjoint Polyakov loop, where C F ¼ ðN 2 c − 1Þ=ð2N c Þ and C A ¼ N c are the Casimir operators of the fundamental and adjoint representations. The parameter δ 8 is the measure of the breaking of Casimir scaling introduced in Ref. [76]. Casimir scaling implies . It was found that δ 8 is small above T c and compatible with zero for T > 300 MeV [76]. Taking into account that we can then write for the subtracted energy In the above equation, we have separated out the factor δV o that breaks Casimir scaling for V o , which is the only purely perturbative input to Eq. (30). This correction was calculated in Ref. [77], Now, using the lattice results for V S and F Q , we can calculate F sub QQ according to Eq. (30) and compare this to the results of direct lattice calculations of F sub QQ . This comparison is shown in Fig. 10 for N τ ¼ 12 at the two temperatures T ¼ 172 and T ¼ 666 MeV. The use of N τ ¼ 12 data for this comparison is justified because they are quite close to the continuum limit. As one can see from the figure, the pNRQCD decomposition formula works very well up to distances rT ≈ 0.3. For T ¼ 172 MeV, the octet contribution is very small and in fact can be omitted. For high temperatures, we could also use F S as a proxy for V s and obtain very similar results. We note that for high temperatures not only the octet contribution but also the breaking of Casimir scaling for V o in Eq. (31) is numerically important. If this breaking would be neglected, then the reconstruction of the Polyakov loop correlator in the right panel of Fig. 10 would only work for rT < 0.2.
We have seen that pNRQCD can explain the lattice results for the Polyakov loop correlator at short distances. Up to distances of about 0.3=T, the behavior of the Polyakov loop correlator can be understood without invoking medium modification of the static quark-antiquark interactions [cf. Eq. (27)]. This is in agreement with our findings for the singlet free energy F S , which also shows that the interaction between the static quark and antiquark is only mildly affected by the medium up to distances rT ≈ 0.3.
In order to compare the free energy F sub QQ with weakcoupling calculations, it is useful to distinguish two regimes. In the regime α s =r ≪ T, we can expand the exponentials in the right-hand side of Eq. (26) and compute them in perturbation theory, obtaining at leading order in α s which is valid up to corrections of order g 5 . Weak-coupling calculations have been performed in this case up to next-tonext-to-next-to-leading order (NNNLO), which is accurate up to order g 7 [31]. In Fig. 11, we compare the NNNLO result for F sub QQ given by Eqs. (54) and (55) of Ref. [31] with the continuum extrapolated lattice result at T ¼ 1600 MeV. We use three values of the renormalization scale: μ ¼ πT, μ ¼ 2πT, and μ ¼ 4πT. The running coupling has been calculated as before. We see that the lattice result and the weak-coupling result agree for rT > 0.2, although the latter has a large uncertainty. At very small distances r, the hierarchy α s =r ≪ T breaks down.
In the regime α s =r ≫ T, the right-hand side of Eq. (26) is dominated by the color singlet contribution. This happens at shorter and shorter distances as the temperature increases. We have seen in Fig. 3 that our lattice data appear to be sensitive to this regime so that we could even define an effective coupling and compare it in Fig. 4 with the strong coupling from the zero temperature force between a static QQ pair.

B. The free energy in the screening regime
In the following, we study to which extent the weakcoupling expression for the free energy is able to describe the lattice data in the electric screening regime. In Sec. IV B, we have shown that the singlet free energy is numerically consistent with the prediction of EQCD at NLO in the range rT ∼ 0.3-0.6. However, studies of the screening regime of the Polyakov loop correlator on the lattice are a daunting task. First, the signal-to-noise ratio of the Polyakov loop correlator is very poor, so we need very large statistics. Second, we find that, due to the smallness of the signal, finite volume effects in the subtracted free energy become important already for rather small separations. Hence, we require lattices with large volumes. We satisfy these constraints by using the new N τ ¼ 4 ensembles with aspect ratio N σ =N τ ¼ 6 and total ensemble sizes of 2 − 6 × 10 4 gauge configurations. The smallness of the signal can be understood from the leading-order expression  of the Polyakov loop correlator. Because of cancellations between singlet and octet contributions, mentioned in the previous subsection, the LO result for F sub QQ in the electric screening regime comes from the exchange of two electrostatic gluons [4,28,29], Therefore, at large distances, F sub QQ is much smaller than F sub S . In Fig. 12, we show our results for the subtracted free energy as a function of rT at different temperatures. The numerical results for different temperatures have been shifted vertically for better visibility. The results at short distances, rT < 0.5, and temperatures below 2 GeV have been extrapolated to the continuum limit (see Appendix C). The continuum extrapolated data are shown as filled symbols in Fig. 12. At large distances, the errors are large for N τ > 4 and thus continuum extrapolations are not possible. To circumvent this problem, we correct the N τ ¼ 4 data for cutoff effects. We define the correction factor K 4 as the ratio between the continuum result for F sub QQ and the Wherever the continuum result is available, the correction factor is known. Reliable continuum results are available only up to distances rT ≈ 0.45. We assume that at distances rT > 0.45 the correction factors are r independent and equal to the correction factor at rT ¼ 0.45. Our analysis shows that cutoff effects are larger at larger distances (see Appendix C). As a consequence, the above procedure gives a continuum estimate that may be higher than the true continuum result (the continuum limit is approached from above). For the highest temperature considered in our study, namely, T ¼ 5814 MeV, no continuum results are available. To make a continuum estimate at this temperature, we note that the correction factor decreases with increasing temperature because the cutoff effects become smaller at high temperatures (see Appendix C). The highest temperature where K 4 can be estimated is T ¼ 1600 MeV. We assume that the cutoff effects at T ¼ 5814 MeV are 2 times smaller than at T ¼ 1600 MeV, i.e., we use Now let us compare the lattice results for F sub QQ with the weak-coupling calculations at r ∼ 1=m D , i.e., with the EQCD calculations. The free energy has been calculated at NLO in EQCD in Ref. [29]. Higher order corrections due to the running of the coupling are numerically important and have to be included in the analysis. These corrections have been computed in Ref. [31]. Therefore, we use Eq. (64) of Ref. [31], together with the NLO result for the Debye mass, Eq. (25), to obtain the weak-coupling result. Keeping only the first term in Eq. (64) of Ref. [31] corresponds to what we call the LO result, evaluated at the scale μ ¼ 4πT. As before, we use ΛM S ¼ 320 MeV and the two-loop running coupling. The weak-coupling results are shown in Fig. 12. The dash-dotted lines correspond to the LO result for μ ¼ 4πT. The NLO results at different temperatures are shown as bands. We vary the scale μ from πT to 4πT to account for the scale uncertainty. The widths of the bands correspond to the scale uncertainties. We see again that the scale dependence reduces significantly at higher temperatures. The NLO result for F sub QQ is larger than the LO result and is significantly closer to the lattice data. Deviations between the lattice data and the NLO result for the scale μ ¼ πT are of the order of the statistical errors. Thus, for distances r ∼ 1=m D , weakcoupling calculations in EQCD give a fair description of the lattice data, at least for the highest temperatures considered. In the next section, we will discuss the very large distance behavior of the static correlators, r ≫ 1=m D , where nonperturbative effects are important.

VI. ASYMPTOTIC SCREENING REGIME
In Secs. IV and V, we have shown that weakcoupling expansions describe lattice results fairly well both at short distances, r ≪ 1=T, and at somewhat larger distances, r ∼ 1=m D , in the electric screening regime, where resummation is necessary. At distances much larger than the Debye length, 1=m D , the magnetic screening, which is inherently nonperturbative, becomes important. Perturbative calculations are not valid in this regime. In this section, we discuss the behavior of the static quark correlators at large distances, i.e., r ≫ 1=m D . We expect that at these distances the correlators will be exponentially screened, i.e., behave like expð−MrÞ=r, with M being a The bands represent the NLO (solid) results, where the scale has been varied as μ ¼ πT, 2πT, and 4πT (solid, dotted, and dashed lines). The dash-dotted lines correspond to the LO result evaluated at μ ¼ 4πT.
nonperturbative screening mass. As we will see, the challenge in the determination of the screening masses is the poor signal-to-noise ratio and possible contaminations from small distance physics.
In Fig. 13, we show our numerical results for F sub QQ and F sub S multiplied by −r as function of rT on a logarithmic scale. At large values of rT, the above quantities should appear as straight lines in the figure. This is indeed the case, i.e., at sufficiently large distances, both F sub QQ and F sub S seem to decay exponentially. In the case of the singlet free energy, the data have some structure: the slope changes as we go from short to large distances. Furthermore, the slope increases with decreasing temperature, implying that the corresponding screening mass in temperature units increases at small temperatures. The fits of the correlators that we use for determining the screening masses clearly show this feature (see Appendix D). The lattice data for −rF sub QQ are relatively featureless. In particular, the slope seems independent of the temperature and the distance range considered. The data for F sub QQ also have much worse signal-to-noise ratio compared to the data for F sub S .
We discuss the details of the analysis of the asymptotic behavior and the extraction of the screening mass in Appendix D. We find that finite volume effects in the asymptotic behavior can be parametrized in terms of a constant shift between the correlators on lattices with different volumes (or aspect ratios N σ =N τ ). After accounting for this constant shift, we find that screening masses are within errors independent of the volume for both F sub QQ and F sub S (see Appendix D). We show the asymptotic screening mass from F sub S in Fig. 14 as a function of the temperature. The screening mass shows some cutoff dependence, which is only mildly temperature dependent. We do not see cutoff effects with a clear dependency on N τ for N τ ≥ 8, although screening masses with N τ ¼ 4 and 6 appear to be systematically smaller. The screening mass is considerably larger than the perturbative Debye mass at NLO [Eq. (25)]. After we rescale m D by a constant A ¼ 1.6-2.0, we find that the temperature dependence of the screening mass and the rescaled Debye mass are very similar for high temperatures, T > 700 MeV. This rescaling factor is consistent with earlier observations in Refs. [70,78] using the LO Debye mass and much lower temperatures. A qualitatively similar enhancement was observed for the electric screening masses calculated from the Landau gauge gluon propagators [74,75]. The two sets of lines in the figure represent the rescaled Debye mass with the scale μ varied from πT to 4πT.
Because of the worse signal-to-noise ratio, the fit of −rF sub QQ is more difficult. We have to restrict the fits to smaller r range. This implies that the contamination from small distance physics could be potentially larger here. On  the other hand, because the slope is independent of rT, we do not see a large sensitivity to the choice of the lower fit range. Our estimate of the corresponding screening mass m QQ is shown in Fig. 15. At high temperatures, this screening mass is clearly larger than the one extracted from F sub S . This is expected and is in agreement with previous findings. We compare our results with the screening mass determination on the lattice from EQCD [79], which is shown as the horizontal band in Fig. 15. Our results and the EQCD result are in good agreement. At temperatures T > 400 MeV, the screening mass from the Polyakov loop correlator is only 25% larger than the naive perturbative result 2m D . This is consistent with the EQCDbased analysis [41]. At temperatures around 200 MeV, we see a decrease in m QQ . Such a decrease in the screening mass has been observed also in the SU(3) gauge theory [17], where it is due to the (weakly) first-order nature of the deconfining transition and thus to the large correlation length of the order parameter (Polyakov loop).
The poor signal-to-noise ratio is clearly a problem for the determination of the screening masses. To overcome this problem, the Wilson lines should be evaluated on smeared gauge configurations before calculating the correlators, as in Ref. [26]. We plan to perform such a calculation in the future.

VII. CONCLUSIONS
In this paper, we have studied the free energy of a static QQ pair and the static meson (singlet) correlator in (2 þ 1)flavor QCD with the aim to analyze at what distances the hot QCD medium modifies the interaction between the static quark and antiquark and whether these modifications could be understood using perturbation theory. We have presented continuum extrapolated results for the free energy of a QQ pair as well as for the singlet free energy in the Coulomb gauge.
At asymptotically high temperatures, gðTÞ ≪ 1, there are three distinct distance regimes: r ≪ 1=T, r ∼ 1=m D ∼ 1=ðgTÞ, and r ∼ 1=ðg 2 TÞ. The physics in these regimes is quite different. In the first case, we expect only small thermal corrections to the quark-antiquark interaction. For r ∼ 1=ðgTÞ, medium effects are significant and the interaction between the static quark and antiquark can be understood in terms of electric screening. At distances r ∼ 1=ðg 2 TÞ, nonperturbative effects due to static magnetic fields are important and the screening is nonperturbative.
Weak-coupling calculations of the Polyakov loop correlator and of the singlet correlator can be organized using effective field theories. For r ≪ 1=T, one can use pNRQCD, while for r ∼ 1=m D , one can use EQCD. At larger distances, the screening masses have to be calculated nonperturbatively using lattice techniques [79][80][81].
For physically interesting temperatures, the relevant coupling is gðTÞ ∼ 1. Thus, the question arises if the above picture and weak-coupling calculations are relevant at all. A related question is how well the different energy scales are separated and at what distances screening sets in. We have addressed these questions through a detailed comparison of lattice results with effective field-theory-based weak-coupling calculations. We have found that the weakcoupling picture provides a fair description for temperatures larger than 300 MeV. For rT ≲ 0.3, medium effects are small, as predicted by pNRQCD. For 0.3 ≲ rT ≲ 0.6, we see significant medium effects and the interaction between the static quark and antiquark is screened as predicted by perturbative EQCD. In this region, screening is controlled by the perturbative parameter m D . At larger separations, nonperturbative effects related to static chromomagnetic fields become important and the weak-coupling approach is no longer valid. In this region, we have determined the asymptotic screening masses by fitting the correlators at large distances.
The above results suggest that weak-coupling pNRQCD could be applicable for studying quarkonium properties at T ≳ 300 MeV. In a broader context, the success of the weak-coupling calculations for the Polyakov loop correlator and the singlet correlator is not completely surprising. Weak-coupling calculations have been found to agree well with lattice results for quark number susceptibilities [59,60,82], equation of state [68,83], and to some extent for the topological susceptibility ( [83,84], see, however, [85]). Furthermore, the observed Casimir scaling of the Polyakov loop expectation values in higher representations [76] is also naturally explained in the weak-coupling picture [73]. The lines correspond to 2 times the NLO Debye mass in temperature units calculated for μ ¼ πT, 2πT, and 4πT (solid, dotted, and dashed) and multiplied by a constant to take into account nonperturbative effects. The open squares correspond to the screening masses determined using N τ ¼ 4 lattices with aspect ratio 6. The horizontal band corresponds to the EQCD result for the screening mass obtained in 3D lattice calculations [79].
As a byproduct of this analysis, we have improved the accuracy of the renormalization constants of the Polyakov loop and the entropy of the single static quark (see Appendix A 3). We have also extended the study of cutoff effects for the T ¼ 0 static energy at short distances (see Appendix B). The latter may be of relevance for improved determinations of the strong coupling α s . . The lattice QCD calculations have been performed using the publicly available MILC code. The data analysis was performed using the R statistical package [86].

APPENDIX A: GAUGE ENSEMBLES
In this Appendix, we give an account of the HISQ gauge ensembles underlying the presented lattice calculations. In Sec. A 1, we discuss the quark mass dependence in the new ensembles and finite volume effects. In Sec. A 2, we cover cuts in molecular dynamics (MD) trajectories due to thermalization, the determination of statistical errors, and the normalization using the MD time histories of the Polyakov loop. We discuss our cuts for the correlators. In Sec. A 3, we finally discuss the renormalization of F QQ and F S .

Lattices
As mentioned in Sec. II, we make use of gauge ensembles from the HotQCD Collaboration and from an earlier publication of the TUMQCD Collaboration (cf. Refs. [7,9,[58][59][60][61]). The parameters of these gauge configurations, the collected statistics that have been used for the correlators in this study, as well as references to the original publications are listed in Tables I-V [87]. This approach avoids small acceptance rates due to large changes of the action during early thermalization and shortens the thermalization phase to typically 20-40 TUs. To be safe, we later omit at least 100 TUs of the initial lattices of each stream in the analysis.
For an unclear reason, we found a problem with thermalization of the N τ ¼ 6, β ¼ 7.825 gauge ensemble. We generated a new ensemble with the same parameters from scratch, which we use in our analysis (cf . Table VIII). In order to take the continuum limit for extremely high temperatures, T ∼ 1-2 GeV, we generated new very fine ensembles with N τ ¼ 12 and 10 (cf . Table VI). In order to ensure control of cutoff effects at the smallest separations, we generated another set of gauge ensembles with N τ ¼ 16 (cf . Table VII). Furthermore, we generated new lattices with a larger light sea quark mass, corresponding to a pion mass of m π ∼ 320 MeV in the continuum limit. These new lattices have    Fig. 16, we show a ratio of the difference between the bare single quark free energies with different light sea quark masses m 1 and m 2 over the combined uncertainty, obtained by adding their statistical errors in quadrature, Δf bare Q fluctuates in a AE1.5σ interval (σ being the standard deviation) and shows no signs of a systematic quark mass dependence. In the middle and right panels, we show screening functions for two representative correlators. The short distance behavior of the screening functions defined as is independent of the quark mass within statistical uncertainties and we do not observe a unique trend towards stronger or weaker screening as a function of the quark masses. Hence, differences can be understood in terms of the much larger statistical errors of the HotQCD ensembles at m l ¼ m s =20, leading to an incomplete cancellation with the squared bare Polyakov loop. From these observations, we have to conclude that the statistical errors on small ensembles may be unreliable. Moreover, we study finite volume effects in the Polyakov loop and in the correlators on a set of ensembles with N τ ¼ 4 and m l ¼ m s =20 or m l ¼ m s =5. We compare gauge ensembles with aspect ratio (AR) N σ =N τ ¼ 4 to gauge ensembles with N σ =N τ ¼ 6 for three different lattice spacings (cf . Table IX). We see small changes in the expectation value of Polyakov loop, which are, given the size of the errors on the N σ =N τ ¼ 4 ensembles, statistically not significant (within one and one and a half statistical errors). We show a comparison of screening functions for different volumes in Fig. 17. The singlet correlators for different volumes are numerically consistent up to rT ≲ 1.25. However, the screening behavior of Polyakov loop correlators for the two volumes differs quite dramatically already at rT ≲ 0.5. Since F sub QQ and F sub S approach in the same volume the same asymptotic constant, the deviation due to different asymptotic values becomes relevant for F sub QQ at smaller separations.

Correlators
As mentioned in Sec. II, we calculate the expectation values and the errors of the correlators with the jackknife method. We generally have one measurement for every ten TUs of molecular dynamics. However, there are some gaps in the MD histories, since some ensembles were not available as coherent MD chains of gauge configurations. We simply ignore these gaps and attach the measurements in ascending order of MD time in each stream. We combine different streams of gauge ensembles into one ensemble for the calculation of the correlators. We omit the first 1000 TUs in the first stream, if the stream did not use a starting lattice with the same lattice spacing. Otherwise, we omit the first 100 TUs in each stream to ensure that correlations between different streams are sufficiently reduced and that the MD histories are not affected by the brief intermezzo of RHMD evolution at the start of a stream. We also combine the Polyakov loop time histories of all streams in the order as the correlators. Since the Polyakov loop is measured after each single TU, there are roughly 10 times as many Polyakov loop measurements as correlator measurements. We divide the Polyakov loop histories evenly over the jackknife bins of the correlators.
We sort data into a fixed number of jackknife bins (30) using bins of equal size and omit leftover measurements. We distribute the Polyakov loop history over an equal number of jackknife bins and divide the correlator on each jackknife bin by the square of the Polyakov loop on the same bin. Since the expectation value of the correlator approaches for asymptotically large r the squared expectation value of the Polyakov loop up to effects due to the finite volume, this ratio, which we call the normalized correlator (see main text) has reduced statistical fluctuations for large separations r in correlators at high temperature. The logarithms of these correlators are the subtracted free energies, which-barring finite volume effects-decrease exponentially for asymptotically large separations r → ∞. We list the maximal separations at which we could obtain a statistically significant signal for F sub S and F sub QQ in Tables I-IX using rT units (usually in the fourth column). For F sub S , we require jF sub S =δF sub S j < 5 for any smaller separation; for F sub QQ , we require at most one instance of jF sub QQ =δF sub QQ j < 2. In some ensembles with insufficient statistics, the cancellation may fail, and thus, we would see a fake signal at large separation. We avoid this by demanding that the modulus of the last accepted points is larger than the modulus of the omitted points and that −F sub as well as screening functions −rF sub and −r 2 F sub have negative slope for sufficiently large rT. We find that this estimator for the largest accessible separations of the QQ pair varies only mildly within about 10% between lattices with the same N τ and similar temperatures as well as similar ensemble sizes. Since the cancellation between the correlator and the squared Polyakov loop is an infinite volume property, incomplete cancellation may as well be a hint of finite volume effects, which can be confirmed for the respective N τ ¼ 4 lattices by inspecting Fig. 17.    If we have data sets with different masses for the same pair of parameters (β, N τ ), we use only the correlator that gives us the best integrated signal-to-noise ratio on the full range that we would use. Hence, our final set of ensembles in the analysis contains two different quark masses.

Renormalization
We determine the renormalization constant 2C Q for the QQ free energies by applying the direct renormalization scheme to the Polyakov loop. This procedure is covered in detail in Ref. [9]. We obtain a set of initial values of C Q from the static energy and include the new T ¼ 0 ensembles with m l =m s ¼ 1=5 [68]. Here, the potential is fixed to the value The Sommer scale r 0 and the scales r 1 and r 2 are defined through dV=drðr 0 Þ ¼ 1.65, dV=drðr 1 Þ ¼ 1.0, and dV=drðr 2 Þ ¼ 1=2. We repeat the direct renormalization procedure of Ref. [9] after including all new ensembles. Because of the availability of higher β at T ¼ 0 and for larger N τ , we do not require iterative direct renormalization anymore. We choose those with smaller uncertainties if multiple ensembles with the same β and N τ are available. In Table X, we summarize the new results for c Q ¼ aC Q with substantially reduced uncertainties. Results for omitted β are commensurate with c Q from Table V in Ref. [9]. As mentioned in Sec. II, we renormalize F QQ and F S by adding twice the renormalized single quark free energy 2F Q to F sub QQ and F sub S . Including the new ensembles, we repeat the extrapolation to the continuum limit (using global fits including N τ ≥ 6 and N −2 τ scaling) (cf. Ref. [9]). In the high temperature region T > 250 MeV, we obtain a continuum limit that is reliable up to T ≳ 2.3 GeV with χ 2 ¼ 97.3 on 95 degrees of freedom after including the full uncertainty of c Q in the global fit. For higher temperatures, the continuum extrapolation is not reliable due to lack of data with N τ > 8. From F Q , we can calculate the entropy of a static quark In Fig. 18, we show the continuum limit of S Q , which agrees with the weak-coupling result at next-to-next-toleading order (NNLO) for T ≳ 1.7 GeV. The values for F Q in the renormalization of the QQ, free energies are taken from the same continuum extrapolation.  . The asymptotic behavior exhibits a strong volume dependence (N τ ¼ 4 and β ¼ 6.664). Cancellation between quark-antiquark free energies and 2F Q is better in the larger volume.  In this Appendix, we discuss the details of the improvement procedure for reducing the cutoff effects at small separations of only a few lattice steps, i.e., r=a < 5, which are due to the explicit breaking of the rotational symmetry in the discretized theory. We develop this procedure for the static energy at zero temperature and apply its result to the quark-antiquark (and singlet) free energies. Lattice results for the static energy are discussed in detail in Refs. [7,58,68]. In Sec. B 1, we analyze cutoff effects in the static energy for the tree-level interaction and the treelevel improvement. In Sec. B 2, we analyze cutoff effects in the static energy for the full interacting theory with the HISQ action and nonperturbative improvement. In Sec. B 3, we apply these corrections to the free energies under study in this paper and give a quantitative account of the nonperturbative improvement that we achieve.

Cutoff effects for the tree-level interaction
Our understanding of the cutoff effects in the static energy of a QQ pair is largely based on the calculation at the tree level, i.e., for one-gluon exchange without a running coupling. We use this to estimate the cutoff effects on the static energy calculated with the HISQ action for β ¼ 10=g 2 ≥ 7.03.
Using the lattice version of the temporal gluon propagator D 00 ðk 0 ; kÞ, it is straightforward to calculate the treelevel static energy on the lattice, which is given by The above equation defines the improved distance r I . Sometimes, in the following, to remark the difference from the bare nonimproved distance r=a, we will indicate the latter with r b =a. In Fig. 19, we show the static energy at tree level calculated for the standard Wilson gauge action and for the improved Lüscher-Weisz (LW) gauge action, normalized by the continuum result of 1=ð4πrÞ. Because of the breaking of rotational symmetry on the lattice, we have to distinguish between on-and off-axis separations. Most of the separations shown in the figure are off axis. For the values of r b =a, which correspond to both on-and offaxis separations (e.g., r b =a ¼ 3), the average of the two is shown. As one can see from the figure, cutoff effects are largely reduced for the improved gauge action and they are very small for r b =a > ffiffi ffi 6 p . We use the tree-level result to reduce cutoff effects in the lattice calculations, namely, for each value of r, we perform the replacement r → r I . This is called the tree-level improvement.

Cutoff effects in the static energy calculated with the HISQ action
The cutoff effects in the static energy with HISQ action have been studied in Ref. [69]. It has been shown that for the HISQ action the cutoff effects follow a similar pattern as in the free theory, namely, cutoff effects are about the same and are very small for r b =a > ffiffi ffi 6 p . Using the tree-level improvement, the cutoff effects can be reduced to 1% or less (cf. Fig. 2 of Ref. [69]). To reduce the cutoff effects even further, we follow the procedure outlined in Ref. [69], but implement new refinements, which are possible due to the fact that lattices for β ¼ 8.0, 8.2, and 8.4 are available by now. We describe these refinements in the following.
In Ref. [58], the static energy has been calculated for lattice spacings down to a ¼ 0.04 fm for m l ¼ m s =20. In  the recent study [68] of the QCD equation of state, the static energy was calculated for lattice spacings down to a ¼ 0.025 fm with m l ¼ m s =5, and it was shown that the quark mass effects on the static energy are negligible at distances r < 0.8r 1 . The lattice spacing was fixed using the scales r 1 or r 2 with r 1 ¼ 0.3106 fm and r 2 ¼ 0.1413 fm. The lattice parameters for the static energy calculations are summarized in Table XI. The static energy calculated on a finer lattice with β ref > β can be used as a reference to estimate the cutoff effects in the static energy on any coarser lattice with smaller values of β. By using the improved distance, we account for the tree-level cutoff effects in this comparison. More precisely, the static energy for r b =a > ffiffi ffi 6 p at β ref provides a continuum estimate of the static energy that lattice data at smaller β can be compared to. To compare the static energies calculated at different β values, one has to use the renormalization constants 2C Q given in Table X. However, the results of the renormalization constants are not sufficiently accurate for comparison of the static energies at subpercent precision. Therefore, we have to shift the reference static energy by a constant within the errors of the renormalization constant for β ref such that the ratio of the potential is close to unity in the distance range r b =a > ffiffi ffi 6 p , where cutoff effects are small. Since we only have the static energy at discrete sets of data points for β ref , we must interpolate between the data using the improved distance r I . The simplest interpolating function is the Cornell potential, VðrÞ ¼ a=r þ σr þ C, and it was used in Ref. [69]. However, now we find that it is not sufficiently flexible to provide accurate interpolations. So we modify the Cornell Ansatz and supplement it by a polynomial in r up to sixth order. We also try a simpler form, namely, VðrÞ ¼ a=ðr lnðkrÞÞ þ σr þ C, which works well in some cases. We also try a spline fit of rV S ðrÞ. We interpolate the data for r b =a ≥ ffiffi ffi 5 p , assuming that cutoff effects are small enough to permit an interpolation with percent level accuracy. In order to extend the available range to smaller r and, hence, to determine estimates of the cutoff effects at smaller separations, we extrapolate towards 10% smaller separations, i.e., down to r b =a ¼ 2.00223. We estimate the continuum limit with the modified Cornell Ansatz and estimate the systematic uncertainty from the difference between fits with this modified Cornell Ansatz and the spline Ansatz. We add the statistical and systematic errors in quadrature.
Some sample results demonstrating this analysis are shown in Fig. 20. We define the correction factors and observe significant deviations of K from zero for r b =a < ffiffi ffi 6 p , although these become smaller for finer lattices. Up to an overall sign (due to normalization) 1 þ K and V lat;free S =V cont;free S are quite similar, although deviations from zero in Fig. 20 are smaller than deviations from one in Fig. 19. In principle, we could divide the lattice data by 1 þ K for any choice of β ref to remove the residual cutoff effects. The finest lattice, i.e., β ref ¼ 8.4 in the present case, can be used as the reference lattice to estimate the cutoff effects in the static energy for any other value of β. In the analysis of Ref. [58], only the finest available lattice with β ref ¼ 7.825 had been used. However, this approach is not optimal for a number of reasons. First of all, the same distance in physical units on a finer lattice corresponds to a larger number of lattice steps than on a coarser lattice. Thus, the signal-to-noise ratio for the same physical distances is worse for finer lattices than for coarser lattices. Second, for larger distances, r ≳ 0.8r 1 , the lattices with smaller light quark mass should be used to avoid systematic errors due to using different quark masses. Finally, the continuum estimates for different β ref suffer from independent statistical fluctuations, which may be reduced by combining different β ref . This leads to a procedure with strong similarities to step-scaling techniques. We vary the reference lattice from β ¼ 8.4 to β ¼ 7.596 and determine estimates for the cutoff effects for lattices with 8.2 ≥ β ≥ 6.488.  Fig. 20.
After determining the β ref -averaged correction factors 1 þ hKi, we put together sets of hKi for different lattice spacings a and interpolate between these. For each r b =a, we use an Ansatz of the form where terms b i0 a 2i can be omitted because the tree-level improvement discussed above removes lattice artifacts of order g 0 a 2 . Moreover, as a term b 0j g 2j 0 would persist in the continuum limit where rotational symmetry is not broken, such a term is forbidden. We perform a bootstrapped fit and vary the number of higher order terms kept in Eq. (B3).
We accept only fits for which we obtain χ 2 =d:o:f: ≤ 2.4 and use their mean as our final interpolation result and the mean of their bootstrap errors as the final statistical error. If we do not obtain fits with reasonable χ 2 =d:o:f., we successively omit the last point at the small β end-the coarsest lattice-and repeat the same set of fits (only for r b =a ¼ 1 we had to omit the last point in order to obtain a smooth interpolation). We estimate the systematic error of the interpolation as the variance of the set of accepted fits. The final χ 2 =d:o:f. is then obtained from the final fit and the input data using an average of the degrees of freedom for the set of accepted fits. We show some sample results of this interpolation in Fig. 21, which describe hKi quite well for all considered lattices. avg.
-0.005 0 0.005 avg. However, we do not always see an unambiguous onset of scaling behavior for very fine lattices. This may indicate that the correction factors hKi do not approach the continuum limit smoothly. For this reason, we conservatively estimate that the last hKiðr b =a; βÞ that we can obtain applies for any larger β and assume a 33% relative systematic uncertainty. For all values of r b =a, we observe a strong increase in the fluctuations of hKi for the coarsest lattices. This effect is due to the choice of the renormalization scheme, in which V S ∼ 0 for such distances. We address this growth of fluctuations by using an alternative renormalization scheme, namely, we shift the static energy by −0.5 GeV with respect to the usual potential scheme. Thus, we defer V S ∼ 0 and the growth of fluctuations to substantially larger distances. In this scheme, we can determine correction factors for the previously problematic range without a large growth of the errors. We repeat the same type of analysis in the shifted renormalization scheme, but use β ref ∈ ½7.030; 7.825 as reference lattices. We observe here that errors seem to be underestimated in the range ða=r 1 Þ 2 ∈ ð0.1; 0.15Þ, which corresponds to β ∈ ð6.5; 6.7Þ. We require that the averaged correction factors hKi have relative errors of at least 5% for β < 6.7 in order to obtain stable fits. We show some sample results in Fig. 22. hKi is significantly smaller in the shifted scheme and the growth of hKi for coarse lattices has been much diminished. In order to compare to the original potential scheme, we have to convert the correction factors back to the original scheme. We show the comparison between both schemes in Fig. 23. The typical differences between both determinations of hKi are about 0.2%-0.4%, but are mostly covered by the statistical errors that remain after converting the shifted scheme back to the standard renormalization. The fair agreement between both schemes indicates that the cutoff effects have been identified appropriately and the systematic error due to scheme dependence is about 0.2%.
If we can determine hKi for given β and r b =a < 4, we assign this systematic error of 0.2% to any data (in the following, denoted as aE), which we add in quadrature to the statistical error. If we fail to obtain hKi at the same β for larger r b =a and at the same r b =a for smaller β, we assume that we have reached the limit of the validity of a controlled determination of hKi. In these cases, we refrain from using the correction factors and instead double the statistical error and add a 1% relative error in quadrature. If the data are too close to zero, i.e., aE < 0.01, we add an additional systematic error estimate of 0.001 absolute to aE. In a case where we cannot determine a correction factor at all (for the coarsest lattices, i.e., β < 6.195), we double the statistical error and add a 2% relative error in quadrature. If the data are too close to zero in these cases, i.e., aE < 0.05, we add the same additional systematic error estimate of 0.001 absolute to aE. At last, we take a weighted average of the corrected results obtained in both schemes using the squared errors as weights. We use the same weights for averaging the errors and add the difference between both as systematical error in quadrature.
Subsuming the aforementioned considerations in our definition of hKi, we define the improved static energy as hKi for fine lattices, i.e., a < 0.075 fm, is usually within AE1%. The first point, r b =a ¼ 1, which also has the strongest tree-level correction, ðr I − r b Þ=r b ≈ −4%, is an exception with hKi ≳ −2%. Combined uncertainties due to statistical errors of the bare static energy and due to systematic errors of the correction factors are usually smaller than 0.5% relative error for fine lattice spacings. Although the latter are the dominant r-dependent source of uncertainty, errors due to renormalization with 2C Q are much larger than those due to correction for the cutoff effects on fine lattices.

Improvement at finite temperature
The improvement procedure that we have developed so far is only relevant for cutoff effects at short distances, where short is meant in terms of only a few units of the lattice step, i.e., typically r b =a ≲ 4. When looking at short distances for finite temperature observables, we usually mean short in units of the inverse temperature, i.e., rT ≪ 1, respectively, r b =a ≪ N τ . In the following, we argue why the same short distance improvement is appropriate for static correlators at finite temperature. At the shortest distances r b =a ≤ 2 cutoff effects usually exceed the small statistical errors of the free energies and cause nonsmooth behavior. This is a major problem for smooth interpolations including short as well as larger distances. In Sec. IV, we have shown that the singlet free energy is vacuumlike up to small corrections for rT < 0.3. Thus, we expect to obtain a similar improvement with the correction factors that have been determined for the static energy at T ¼ 0. Namely, we define However, for N τ ≤ 6, it is not obvious if these improvement procedures are appropriate, since r b =a ∼ 2-3 is already within the electric screening regime. We show a few sample singlet screening functions with tree-level or nonperturbative improvement and without improvement in Fig. 24. For N τ < 8, modifications due to improvement are small compared to medium effects and do not alter the result significantly. On the finer lattices, i.e., N τ ≥ 8, improvement smoothens the short distance regime for r b =a < ffiffi ffi 6 p and substantially helps with obtaining interpolations of the data. The systematical uncertainties of the improvement factors are usually at least as large as the statistical errors of the finite temperature results at short distances. These larger errors are beneficial for smooth interpolations of the improved data. Since we use the same tree-level and nonperturbative corrections for V S and F S , we expect that leading corrections are subject to strong cancellations in the difference V S − F S . We simply use unimproved distances r b =a and omit nonperturbative correction factors.
In Sec. V, we have shown that the Polyakov loop correlator can be understood through its decomposition into color singlet and octet in pNRQCD for rT < 0.3 at all temperatures in the study. Using Eq. (30), we can obtain a correction prescription for F QQ on the basis of the correction factors hKi for the T ¼ 0 static energy. We write the improved normalized Polyakov loop correlator C sub P;I in terms of the improved static energies, V S;I neglecting any Casimir scaling violations, We know V S;I ¼ V S =ð1 þ hKiÞ and we can use hKi ≪ 1 to expand C sub P;I in hKi, eventually obtaining The need for corrections and the validity of the pNRQCD formula are both restricted to short distances. In this range, we have seen that the reconstruction using F S works as well as using V S , but provides access to the full β range. Hence, we replace V S by F S to obtain our final expression which is valid to leading order in hKi. We show a few sample screening functions with tree-level or nonperturbative improvement and without improvement in Fig. 25. The effects of improvement are hardly relevant for F QQ on most lattices. In the cases where these effects are visible, the corrections generally lead to a smoother free energy at short distances.

APPENDIX C: INTERPOLATIONS AND EXTRAPOLATIONS
In this Appendix, we discuss the details of the interpolations of the subtracted quark-antiquark (and singlet) free energies and how we extrapolate F sub QQ (and F sub S ) to the continuum limit. We do not repeat here the discussion of the physical underpinning of the presented analysis, which we discuss in detail in Sec. III A. Here we also discuss the interpolations and extrapolations of In the following, we refer to either F QQ or F S just as the free energy F unless we apply different treatment to both and use these quantities in temperature units, namely, we fit f sub ¼ F sub =T and ½V S − F S =T. This Appendix also covers the calculation of the r derivative of F. We usually propagate statistical errors through the fits by generating 100 bootstrap samples of mock data assuming Gaussian errors and perform the fits on the mock data. We perform fits using either linear regression, the built-in smoothing splines of the R statistical package [86], or use for nonlinear fits the nlme package [88].
We arrive at our continuum result through the following procedure. (1) We consider the subtracted free energies F sub defined in Eqs. (10) and (11) and V S − F S . We omit data at large rT due to signal loss (cf. Appendix A) and perform the tree-level and nonperturbative improvement if needed (cf. Appendix B). (2) We interpolate in rT and T using two different datadriven schemes. We use local fits in the first scheme and global fits in the second scheme. (3) We calculate weighted averages of the results from the two schemes and estimate systematic errors due to model dependence. (4) We extrapolate the combined results to the continuum limit for each individual ðrT; TÞ point using different subsets of the N τ and different Ansätze for scaling behavior. (5) After inspecting scaling behavior for sets of representative ðrT; TÞ points, we compile the final continuum result using different extrapolations in those regimes where they are appropriate. (6) Finally, we obtain the renormalized continuum result by adding 2F Q (cf. Appendix A). As stated, we follow a data-driven scheme in fitting the rT and T dependence. Namely, we iteratively try to find the most simple, smooth function of rT that incorporates the main features known from weak-coupling results (LO short distance behavior and exponential damping) and describes the largest subset of the statistically relevant data. We reduce the fit range at the large distance end in small steps if we fail to fit the data with any reasonable choice of smooth function. If the data set becomes too short, the statistical errors of the data at the shortest distances and of the nonperturbative corrections must have been underestimated. This implies that data are nonsmooth at the shortest distances. In these cases, we have to enlarge the errors of all data by a small amount that permits a smoother curve and restart the procedure. We end the procedure when meeting a stopping criterion in terms of both the residues at the shortest distance and globally and in terms of a reasonable determination of the propagated errors of the interpolating function. We discuss the stopping criteria for each of the fits in the following sections.

Local fits
We perform local interpolations first in the quarkantiquark separation and later in the temperature. For the former, we explicitly account for correlations between data at different separations obtained within the same gauge ensemble by fitting on the individual jackknife bins. We have tested that a different ordering of the fits, i.e., first interpolating in the temperature for fixed distances rT and then interpolating in the distance rT for fixed temperature, does not change the results in a statistically significant way. However, the latter ordering does not permit one to account for the correlations within each ensemble.

a. Distance interpolations
We independently interpolate F sub ðr; T; aÞ or rather F sub ðrT; β; N τ Þ as a function of rT for each N τ and each β. We split the data into intervals of small and large separations, where we treat the first 16 points as data at small separations and the rest including some overlap region as data at large separations. Before interpolating data at large separations, we calculate residues at large separations by extrapolating with the fit to data at small separations. The residues are computed with data at large separations that did not enter into the fit at small separations. If these residues satisfy χ 2 =d:o:f: < 1 (d:o:f. being the data at large separation), we forgo the separate interpolation at large separations and just use the extrapolation. If we have separate interpolations for small and large separations, we match the two interpolations smoothly.
For the local fits, we assume smoothness and screening of the subtracted free energies and use a modified, screened Cornell potential (MSCP) Ansatz We keep a parameter C to account for a possible offset due to incomplete cancellation in the subtraction of 2F Q , but also try the same fits with C ¼ 0, i.e., omitting the parameter altogether. We demand that the screening mass parameter M is positive, but also try fits with M ¼ 0 for low temperatures, i.e., T < 235 MeV. We generally use n min ¼ 0 and limit the number of monomials by n max ≤ 8 if enough data points are available. If we fit for small separations with at least six data points and not very high temperatures, namely, with at most T ≤ 500 MeV, we limit the minimal number of monomials as n max ≥ 1. Otherwise, we use n max ≥ 0. In the denominator, we choose at short distances α ¼ 2 for F sub QQ due to the form of the LO result [Eq. (32)], and α ¼ 1 for F sub S . At larger separations, we use α ¼ 1 both for F sub S and F sub QQ . Lastly, we use α ¼ 0 for V S − F S . In the latter case, we also limit the range to rT ≤0.5 for N τ > 4 and to rT ≤ 0.75 for N τ ¼ 4. Because of the small size of the gauge ensembles for N τ ¼ 16, we enlarge the corresponding statistical errors by a 1% relative error estimate.
We vary n max and thus perform nonlinear fits for different polynomial orders N n ¼ n max − n min þ 1 in order to choose the most simple Ansatz that matches the data for short or large separations. In a first cycle, we interpolate only the jackknife average to determine an appropriate fit function that permits meeting a stopping criterion. In a second cycle, we perform the interpolation on all jackknife bins, thus taking correlated fluctuations of data into account and propagating the statistical errors. If the need arises, we readjust the fit function as in the first cycle. We estimate the bias from the difference between the average of the fits on the individual bins of data and the fit on the average of the bins of data and add the bias to the statistical error in quadrature. As stopping criterion of the loop, we demand that the fit converges on at least 20% of all bins with χ 2 =d:o:f: ≤ 2.4 5 and that-after we exclude fits on bins with χ 2 =d:o:f: > 2.4 from an average-the average of the fits satisfies χ 2 =d:o:f: ≤ 1.2. If we reach small enough χ 2 =d:o:f. for some n max , we also attempt a fit with polynomial of order n max þ 1 with the constant term C omitted. Between these, we choose the fit with smaller χ 2 =d:o:f.
We perform a nested loop over different fits starting from the minimal permitted n max and increase towards the maximal permitted n max in an inner loop. If we fail to meet the stopping criterion, we omit points and restart the inner loop. For small (large) separations, we omit points at the large (small) r end of the fit range. If we fall below nine data points in this process, we restore the full fit range, but add a systematical error estimate on top of the previous errors to the data. We add 0.5% times the minimal χ 2 =d:o:f. of the recent inner loop as relative error estimates and restart the inner loop. If we fail to fit for small (large) separations with at most 20% (30%) relative error added, we consider the fit as a failure. We restart the outer loop without a screening mass parameter M for small separations at low temperatures, i.e., T < 235 MeV, reasoning that we do not have enough data to resolve screening as well as subleading terms in the polynomial. For small separations, we aim at fits with sufficiently small residue on the first point. If we cannot achieve χ 2 ðr=a ¼ 1Þ < 1, we restart the outer loop and permit χ 2 ðr=a ¼ 1Þ < 1.5 or even 2. We can always get acceptable fits under the latter condition. We use the fits for large separations to extrapolate up to the first point, where the fit to F sub becomes compatible with zero within errors. We test for model independence of the fit procedure by performing smooth spline fits to rF sub ðrTÞ fits with 100 bins of bootstrap mock data treating all errors as independent. For the latter, we use as few knots as possible and a smoothing parameter as large as permissible to still obtain χ 2 =d:o:f: ≤ 1.2. We find fair agreement between all these fits within their respective uncertainties, albeit with underestimated errors for the bootstrapped nonlinear fit and with overestimated errors for the spline fit. We consider the MSCP fits on the jackknife bins as our final local rT interpolation of F sub . We obtain the rT derivative analytically for every fit function and insert the fit parameters on the individual bins to compute the average and propagate the errors.
For V S − F S , we follow more or less the same procedure. We omit the screening mass altogether and do not use jackknife bins, since we obtain V S − F S by subtracting jackknife averages of statistically independent data sets whose sizes may differ substantially. In practice, the actual fluctuations appear to exceed the statistical errors, which have been added in quadrature. Since many data for V S − F S are very close to zero, we demand that absolute errors 6 are at least 0.001T. In the case of enlarging the errors, we add 0.002T times the minimal χ 2 =d:o:f. of the recent inner loop as absolute error estimates and restart the inner loop. We do not extrapolate to larger separations here. Other than these, the local fits to V S − F S follow the same procedure as the fits to F sub .

b. Temperature interpolations
In the next step, we independently interpolate the MSCPinterpolated results for f sub ¼ F sub =T or ½V S − F S =T in T for each rT and fixed N τ . We perform separate temperature interpolations for the rT derivatives. We use smoothing spline interpolations that penalize higher derivatives. We require for the internal smoothing parameter (SM) of the built-in smooth splines of the R statistical package [86] SM ≥ 0 for f sub and SM ≥ 0.2 for the derivatives. We limit the maximal smoothing parameter SM for the temperature interpolations to 0.45 and include data for all available temperatures. We perform separate interpolations in overlapping low (T < 255 MeV) and high (T > 180 MeV) temperature intervals. If the fits yield χ 2 =d:o:f: > 1.2, we enlarge the errors of the input by a factor ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi χ 2 =d:o:f p . Here we have to stress that the input for these interpolations are intermediate mock data, whose errors may have been underestimated. If the fits in these intervals fail to produce χ 2 =d:o:f: < 3.6 7 without enlarging the errors, we remove the highest (or lowest) temperature in the overlap region from the fit and repeat the fits. We always keep at least one temperature common between both fits. This procedure is necessary since data at very low temperatures usually have much larger uncertainties than at higher temperatures. Thus, a global interpolation cannot resolve the increase of f sub towards low temperatures properly. In agreement with predictions of weak-coupling results (cf. Secs. IV and V), the T dependence is much flatter at higher temperatures than at lower temperatures. In the overlap region, the second derivative is quite large, and thus, use of a smoothing parameter penalizing higher derivatives is hardly appropriate. In the full overlap region, we match the fits for low and high temperatures smoothly.
Since V S − F S ∼ 0.02T for very small distances and high temperatures, see Sec. IVA, this enlarged error cannot create fake plateaux at short distances. 7 This restriction is heuristic. When enlarging the errors after obtaining χ 2 =d:o:f: > 3.6, the fits generally fail to reproduce the input values in the overlap region well.

Global fits
We perform separate global fits to the free energies and to V S − F S for each N τ . For the free energies, we perform the fits over the full temperature range as well as for low or high temperature intervals with T ≤ 235 MeV or T ≥ 175 MeV, respectively. Because of the small size of the gauge ensembles for N τ ¼ 16, we enlarge the corresponding statistical errors by a 1% relative error estimate. For V S − F S , we limit the rT range to rT ≤ 0.5 for N τ > 4 and to rT ≤ 0.75 for N τ ¼ 4, but fit all available temperatures at once. In order to obtain stable fits, we have to enlarge the errors of V S − F S to 5% relative error or 0.001T absolute error, whichever is larger.
For the global fits, we use a generalized temperaturedependent modified screened Cornell potential (GTDMSCP) Ansatz. In the GTDMSCP Ansatz, we generalize Eq. (C1) as We did not observe any systematic temperature dependence in the values of C from local fits using the MSCP Ansatz [Eq. (C1)]. Hence, we must omit C in the global fits and consider effects leading to C ≠ 0 in the local fits as temperature-dependent uncertainties that have to be averaged out in the global fit. In the GTDMSCP fits, we assume that the dominant temperature dependence is in the double polynomial, while we assume that the temperature dependence of the screening mass MðTÞ is rather mild. 8 We model the temperature dependence of the screening mass as where we restrict the range of M 1 to the interval ð−1; þ1Þ and normalize the temperature logarithm by T 1 ¼ 1 GeV to ensure a relatively mild temperature dependence of MðTÞ. This choice of T 1 restricts the logarithm to about ð−2; þ2Þ, which is sufficiently flexible to permit even negative screening masses at one end of the temperature range. As we typically obtain −0.5 ≲ M 1 ≲ 0.5 in fits, our model is sufficiently flexible. Unless it comes out as zero, M 0 is generally positive and larger than jM 1 j by an order of magnitude.
To model the dominant temperature dependence in the GTDMSCP Ansatz in a very general form, we use sequences of consecutive powers of the temperature T m i and vary the lowest and highest powers. We use m i between m min ≥ −4 and m max ≤ þ3 and sequences with lengths N m ¼ m max − m min þ 1 ¼ 1; 2; …; 6. After rescaling the temperatures as T=T 1 , we obtain mostly coefficients of natural size. We construct all consecutive sequences including m i ¼ 0 and iteratively increase N m when needed. We use n max ¼ 8 as the highest power of rT, but observe that the global fits generally favor lower powers of rT. As leading inverse power of rT, we use α ¼ 1 for F sub S and α ¼ 2 for F sub QQ . The underlying reasoning for this choice is the form of the LO result [Eq. (32)].
We perform nonlinear fits for different choices of N m and N n to determine the most simple Ansatz that matches the data. We use three nested loops for searching the optimal fit function. The two outer loops are the same as for the local fits, while the innermost loop is the loop over the sequences of consecutive powers of T. In a first cycle, we perform the global fits on the jackknife average to identify an appropriate fit function. In a second cycle, we use bootstrap mock data to propagate the statistical errors. If the need arises, we readjust the fit function in the second cycle again. If we fail to reach convergence, we add 0.5% times the minimal χ 2 =d:o:f. of the recent loop and restart the outermost loop. We observe that we almost universally need sequences from T −3 to T 0 and sometimes in addition T 1 . For this reason, we eventually start from sequences of length m ¼ 4 for our final fits. With regard to the rT dependence, we demand for the free energies at least two and at most three terms beyond the leading powers. We summarize the global fits of F sub S and F sub QQ in Tables XII and XIII. Requirements for powers of the rT polynomials, enlarged errors, and achieved χ 2 =d:o:f. results are quite similar for the local fits.
For V S − F S , we omit the screening mass completely and demand at least three terms beyond the leading power, which we fix through α ¼ 0. As discussed in the subsection on local fits, V S − F S appears to be affected by fluctuations that are larger than the statistical errors. In order to obtain stable interpolations, we have to enlarge the errors globally to the larger of 0.001T absolute error or 5% relative error before fitting. We summarize the global fits of V S − F S in Table XIV.    In Fig. 26, we show a comparison of the low and high temperature global fits and the corresponding local fits before temperature interpolation with the underlying data for N τ ¼ 12. Within the range of fitted data for T ≤ 220 MeV we obtain fair agreement between local and low temperature global fits, but not with the other global fits. For T ≥ 220 MeV, we find fair agreement between local and global fits for high temperatures or on the full temperature range. Therefore, we switch for T > 220 MeV from the low temperature fits to the high temperature fits. We estimate systematical errors of the global fit from half of the difference in the range T ∈ ð190; 230Þ MeV and add it to the statistical error. We obtain the rT and T derivatives analytically for every fit function and insert the fit parameters on the individual bins to compute the average and propagate the errors.

Average of local and global fits
The assumption of smoothness in the temperature leads to substantially smaller errors in the interpolating functions than in the actual data (see Fig. 26). Hence, we assume that the errors of the interpolating function may be underestimated and we have to account for that. This is particularly relevant at the shortest and largest distances. At the shortest distances, there are only a few points with relatively large systematic errors, which have only little weight in the global fits. At the largest distances, the local fits are oblivious of the smoothness condition between neighboring temperatures and are sensitive to incomplete subtraction. These complementary systematic uncertainties can be directly observed in Fig. 26.
We eventually calculate the weighted average between the local and global fits before the continuum extrapolation. The former do not require that the subtracted free energies approach zero at large rT, while the latter do not require that the screening mass parameter is positive. Hence, they have different systematic uncertainties for large rT, which can be controlled better through an averaging procedure. We estimate the systematical error from the difference between local and global fits. We restrict the growth of the error of the averaged result such that the 1σ band of the averaged result is contained by the combined 1σ bands of the local and global fits. Moreover, to account for any  remaining uncertainties at short distances, we demand that the relative error must not be smaller than the difference between relative errors of the corrected and uncorrected data for the nearest point in the ðrT; TÞ plane. This is still not sufficient for a reliable continuum extrapolation. We find that quite often interpolated data with N τ ¼ 10, 8, and sometimes 6 are in the scaling window, while some interpolated data for N τ ¼ 12 or 16 are not consistent with this scaling behavior (see the scaling plots in Sec. C 4). This is direct evidence that interpolations with N τ > 10 still have underestimated uncertainties. For this reason, we enlarge the errors of interpolations by hand for N τ > 10. We add to interpolations with N τ ¼ 16 a 1% relative error for rT < 0.2 and a 2% relative error for 0.2 ≤ rT < 0.3. We add to interpolations with N τ ¼ 12 a 0.5% relative error for rT < 0.2 and a 1.5% relative error for 0.2 ≤ rT < 0.3.

Continuum extrapolations and error estimates
We extrapolate the subtracted free energies F sub QQ and F sub S , the derivative ∂F S =∂ðrTÞ (and thereby α QQ ½F S ), and V S − F S to the continuum limit. With all continuum extrapolations, we attempt to cover as much of the ðrT; TÞ plane as we can control with reasonable uncertainties. In particular, since we are interested in the different regimes of the QQ interactions, we have to anticipate a drastic variation of the cutoff effects of each quantity depending on the region of the ðrT; TÞ plane. Moreover, depending on the region, we have very different subsets of the data available. Namely, for small separations, we can use only fine lattices with large N τ , and  for large separations, we have only data with smaller N τ due to the increasing noise problem with larger N τ . For these reasons, we must assume that a global continuum extrapolation may fail to reflect the associated uncertainties of the data as well as the richness of the underlying physical structures appropriately. Instead, we perform independent continuum extrapolations using the averaged, interpolated results for each point of the ðrT; TÞ plane.
Because of the structure of the HISQ/tree action, the leading cutoff effects are Oðα s a 2 ; a 4 Þ. Since we take the continuum limit at fixed temperature and T ¼ 1=ðaN τ Þ, this translates into a parametrization of leading cutoff effects through inverse powers of N 2 τ . We use the Ansatz and vary the set of N τ values that are included in the fit. Furthermore, we either use just the a term (assuming data are already in the continuum limit within errors) or include the b=N 2 τ term, assuming a 2 scaling. In extrapolations that include data with N τ ¼ 6 or 4, we also include the c=N 4 τ term. As we perform independent continuum extrapolations, we do not make any assumptions about the ðrT; TÞ dependence of a, b, and c. For this reason, our continuum extrapolations do not necessarily yield a smooth function in the ðrT; TÞ plane. We list the different types of fits used for continuum extrapolation in Table XV. If the number of data with different N τ exceeds the number of parameters of the interpolating function, we perform a linear regression fit and directly propagate the errors. After propagating the errors through the linear regression, we rescale them by 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi χ 2 =d:o:f p . in cases where we obtain χ 2 =d:o:f: < 1. Occasionally, this may lead to large spikes in the error of the continuum result if the underlying data are aligned particularly well due to fluctuations. Otherwise, we estimate the errors of the continuum limit with bootstrap mock data. We instead omit c=N 4 τ if we have only three different N τ available.
We have to consider the possibility that the sequential interpolations or the global fits together with the average and systematic error estimates may still underestimate the actual uncertainties of the correlators. If necessary, we enlarge the errors (of the averaged data underlying the continuum extrapolation) to reach χ 2 =d:o:f: ≤ 1.2. For each datum, we determine in the ðrT; TÞ plane the closest point of the nonperturbatively corrected data set for the same N τ . We permit enlarging the relative error of the averaged result up to the relative error of this closest corrected original datum. We simultaneously enlarge errors of the results for all N τ in four steps up to the permissible maximum for each and repeat the fit. In Figs. 27-30, we show the approach to the continuum limit for the various free energy observables, usually in terms of a ratio of the lattice result over the continuum limit. We show two representative distances: one in the vacuumlike regime and one in the electric screening regime. We have shifted the ratios for different temperatures vertically for better visibility. The colored horizontal lines indicate the continuum value for each temperature to facilitate reading off the size of the cutoff effects for each N τ .
We show the approach of F sub QQ to the continuum limit for a few representative pairs ðrT; TÞ in Fig. 27. In general, we see that within an a 2 -scaling window F sub QQ decreases towards the continuum limit and that a 4 -scaling contributions come with opposite sign. Within numerical uncertainties, F sub QQ with N τ ¼ 6 is always inside the a 2 -scaling window for T > 180 MeV, while F sub QQ with N τ ¼ 4 is never inside the a 2 -scaling window. In general, we see that the approach to the continuum limit becomes steeper for lower temperatures and larger separations and flatter for higher temperatures and shorter separations, i.e., cutoff effects become larger from the upper left to the lower right corner of the ðrT; TÞ plane. For small distances up to rT ∼ 0.15, cutoff effects are at most a few percent or not even unambiguously resolved at all. For intermediate distances up to rT ∼ 0.30, cutoff effects grow up to about 20% in the transition region and 5% for temperatures beyond 1 GeV. For large distances beyond rT ∼ 0.30, cutoff effects are even larger, 30% and more for N τ ¼ 6, and stay at such level even for the highest temperatures. We summarize in Table XVI which continuum extrapolations we combined in our final result for F sub QQ .
We show the approach of F sub S to the continuum limit for the same representative pairs ðrT; TÞ in  with N τ ¼ 6 is always inside the a 2scaling window for T > 180 MeV, while F sub S with N τ ¼ 4 is never inside the a 2 -scaling window. In general, we see as for F sub QQ that the approach to the continuum limit becomes steeper for lower temperatures and larger separations and flatter for higher temperatures and shorter separations, i.e., cutoff effects become larger from the upper left to the lower right corner of the ðrT; TÞ plane. For small distances up to rT ∼ 0.15, cutoff effects are at most a few percent or not even unambiguously resolved at all. In particular, the coefficient b seems to switch to the opposite sign above T ∼ 600 MeV and rT ≪ 0.2. We consider this as an indication that these data may be considered in the continuum limit within errors. For intermediate distances up to rT ∼ 0.30, cutoff effects grow up to about 15% in the transition region and 4% for temperatures beyond 1 GeV. For larger distances beyond rT ∼ 0.30, cutoff effects are even larger and gradually increase up to 40% for N τ ¼ 6 at rT ∼ 0.60. For the highest temperatures cutoff effects are about half this size. We summarize in Table XVII which continuum extrapolations we combined in our final result for F sub S . We show the approach of α QQ ½F S to the continuum limit for the same representative pairs ðrT; TÞ in Fig. 29. Contrary to f S we see that within an a 2 -scaling window α QQ ½F S increases towards the continuum limit at short distances, rT ≪ 0.4, while it tends to decrease towards the continuum limits for larger distances. This change in behavior is clearly due to the physical regimes, namely, the vacuumlike regime and screening regime. The a 4scaling contributions always lead to an increase towards the continuum limit. Within numerical uncertainties, α QQ ½F S with N τ ¼ 6 is always inside the a 2 -scaling window. Even α QQ ½F S with N τ ¼ 4 is close to the a 2 -scaling window for the smallest separations. In general, we see for α QQ ½F S that the approach to the continuum limit is different in the vacuumlike regime and in the screening regime. In the vacuumlike regime, it seems to become slightly steeper for higher temperatures and some intermediate separations, while in the screening regime, it behaves rather in the opposite way. For small and intermediate distances up to rT ∼ 0.30, cutoff effects are at most about 5% percent or not even unambiguously resolved at all. For larger distances beyond rT ∼ 0.30, cutoff effects eventually change their sign and become even larger and similar to those of f S in terms of sign and magnitude. We summarize in Table XVIII which continuum extrapolations we combined in our final result for α QQ ½F S .
Lastly, we show the approach of V S − F S to the continuum limit for a similar set 9 of representative pairs ðrT; TÞ in Fig. 30. We cannot plot the ratio of results at finite lattice spacing and the continuum result, since data are actually crossing zero at different points in the ðrT; TÞ plane. In general, we see that within an a 2 -scaling window V S − F S decreases towards the continuum limit and that a 4scaling contributions come with the opposite sign. Contrary to F sub S , we see that the a 2 -scaling window extends only to N τ ¼ 8. Although V S − F S with N τ ¼ 6 appears to be systematically outside the a 2 -scaling window, the a 4scaling contribution is rather mild. In general, we see for V S − F S that the approach to the continuum limit appears very flat for low temperatures, T < 200 MeV, and quite temperature independent for T > 200 MeV. Cutoff effects become systematically larger for larger separations. For small distances up to rT ∼ 0.15, cutoff effects are not even unambiguously resolved and data at finite N τ fluctuate around the continuum result without clear scaling behavior. We consider this as an indication that these data may be considered in the continuum limit within errors. For intermediate distances up to rT ∼ 0.20 and temperatures above 200 MeV, cutoff effects grow in size from 0.01T to 0.03T with increasing separation and show clear scaling behavior. For larger distances around rT ∼ 0.30 and beyond, cutoff effects seem to stay at a similar level. We summarize in Table XIX which continuum extrapolations we combined in our final result for V S − F S .
We universally find that cutoff effects at the shortest distances that we can access with multiple N τ are consistent with zero within errors. Thus, we extend the continuum TABLE XVIII. Final continuum extrapolation of α QQ ½F S . We list rT and T ranges in the first column, suitable extrapolations in the second column, and weights for the difference (if applicable) as systematic error estimates in the third.
rT; TðMeVÞ ½ minðN τ Þ; maxðN τ Þ; p Δ syst result to even shorter distances by taking the results for the (respective) largest N τ in the extremely short range as a continuum estimate. These continuum estimates are included in the continuum results that are plotted in the figures in the main body of the paper.

Continuum estimates for T = 0 results
For the continuum estimates of the zero temperature static energy, we use the nonperturbatively improved results at finite lattice spacing. We use the lattices in Table XI and the improvement procedure discussed in Appendix B.
After nonperturbative improvement, the results at different lattice spacings are consistent within their respective errors. Therefore, in the range where we have results available at multiple β, each of these results is within its uncertainties in the continuum limit. For the shortest distances, r < 0.029 fm, we assume that the nonperturbatively improved results on the finest lattice are consistent with the continuum limit and treat these as our continuum estimate. Instead of following a similar approach as we did for T > 0, (cf. Sec. C 4), we construct our entire continuum estimate for T ¼ 0 from nonperturbatively improved results at finite lattice spacing. The three reasons for this change of procedure are that we do not have to consider T-or rTdependent cutoff effects, that the results on finer lattices anyways do not contribute strongly to the continuum limit at larger distances, and that the finest lattices have a larger light sea quark mass than the coarser lattices. The second reason can be understood from the fact that such distances correspond to rather large values of r=a on finer lattices, and thus, results on finer lattices have accumulated more noise at the same distance. We construct our continuum estimate by using only points with r=a ≲ 6 for each lattice spacing to keep the errors small. Quark mass dependence as the third reason has been discussed already in Appendix B. By restricting the lattice results with larger quark mass to distances much smaller than 0.8r 1 , we avoid that the quark mass dependence is numerically important.
For the T ¼ 0 results in the figures, we obtain two bands as follows. We use results with five different lattice spacings, namely, with β ¼ 8.4, 8.0, 7.825, 7.28, and 6.74. We plot the nonperturbatively corrected results for V S for each lattice spacing in a different r range to ensure that we always use data at only a small number of lattice steps, r=a ≲ 6. For lattices with m l ¼ m s =5, we use even shorter separations, r=a ≲ 4. For the coarser lattices, we leave out the first few points, where cutoff effects are somewhat large, and thus, the nonperturbative improvement introduces a larger systematic error. The reason for these short distance cuts is that we aim at obtaining a smooth error band. We summarize this construction in Table XX and show the corresponding T ¼ 0 result with black pentagons and a black band in Fig. 1. For the T ¼ 0 result in Figs. 2 and 4, we use a derivative of this T ¼ 0 result obtained numerically from a fit with an Ansatz as a modified Cornell potential (see Appendix B for a discussion of these fits).
However, this is still not enough for a comparison of F QQ from high temperatures to the T ¼ 0 limit. For this comparison, we make use of information obtained from the direct comparison between V S and F S in Sec. IVA. Figure 6 shows that ½V S − F S =T is for sufficiently small rT within numerical uncertainties constant and independent of the temperature for T ≳ 300 MeV. This constant is small compared to the numerical uncertainty introduced due to the renormalization constant 2C Q at high β values. Hence, within the numerical errors, the renormalized result for F S is a fair estimate of V S as long as T ≫ 300 MeV and r ≪ 0.2=T. We employ the continuum result for F S in the range 0.09 ≤ rT < 0.17 for T ≥ 800 MeV, namely, for the temperatures T ¼ 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, and 2.2 GeV, as an estimate for the zero temperature result at even shorter distances, reaching down to r ∼ 0.01 fm. This estimate for V S is shown as a light gray band in Fig. 1. In this Appendix, we discuss the details of the analysis of the asymptotic behavior of the subtracted free energies and the extraction of the asymptotic screening masses. For this analysis, we perform fits on all jackknife bins to account for correlated fluctuations of neighboring data points. The discussion of finite volume effects in Appendix A revealed that screening functions for static quark-antiquark free energies are quite susceptible to finite volume corrections. Thus, we fit F sub for large separations as assuming that the one-particle exchange is the dominant process for large separations. The asymptotic cancellation between the squared Polyakov loop and the static quarkantiquark correlators strictly holds only in infinite volume. Hence, we account with the fit parameter C for possible offsets due to an incomplete cancellation with 2F Q at finite volume. We determine C by fitting a constant to the asymptotic tail of F sub . We define this tail as made of all points with signal-to-noise ratio less than one and of all points at larger rT. By subtracting the parameter C from F sub , we account for the asymptotic finite volume corrections and bring results with different volumes into much better numerical agreement. For β ¼ 6.664, we show the effect of the correction on F sub QQ in Fig. 31. There, we obtain a good signal in the larger volume up to rT ≈ 1.3. After the finite volume correction, the central value of the result in the smaller volume sits on the top of the finite volume corrected result in the larger volume (albeit with large errors that render analysis of the last few data points moot). For different β, the picture is quite similar. With the singlet free energy, we can access much larger ranges, but the main conclusions stay the same (cf. Table IX for the uncorrected results).
Thus, we conclude that we may consider the other two fit parameters A and M as volume independent within errors. We determine ln½rTðF sub − CÞ ¼ − ln½AT − MrT and fit in the range up to the first (uncorrected point with signal-tonoise ratio less than one) through linear regression. By varying the minimal and maximal separations ðrTÞ min and ðrTÞ max that are included in the fit, we obtain a local definition of the two parameters A and M on each jackknife bin. We propagate the errors on each bin via regression and enlarge the errors by 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi χ 2 =d:o:f p . if χ 2 =d:o:f: < 1. We neglect results where regression errors of either parameter are 100%. We propagate the statistical errors with the jackknife method and add the estimate of the bias in quadrature. Finally, we calculate for each ðrTÞ min a weighted average of the parameters obtained with different ðrTÞ max and estimate the error by adding the average jackknife error and an estimate for the bias due to dependence on ðrTÞ max in quadrature. We interpret the fit parameter M as the screening mass in units of the temperature.
We show the local screening mass m QQ =T from F sub QQ in Fig. 32. There is fair numerical agreement between the screening masses for both aspect ratios. For the smaller  volume, we fail to determine m QQ =T reliably for the highest temperature, which is caused by the much smaller signal and, thus, the much worse signal-to-noise ratio. At lower temperatures, the agreement between m QQ =T for different volumes is quite remarkable, in particular, the screening mass assumes a plateau already at ðrTÞ min ∼ 0.4 − 0.5, which is approached from above. We observe that the difference between m QQ =T for ðrTÞ min ¼ 0.25 and for ðrTÞ min ¼ 0.5 is about 0.5. We use this difference to correct for screening masses m QQ =T determined at ðrTÞ min ¼ 0.25 on lattices with aspect ratio N σ =N τ ¼ 4 and assign a systematic uncertainty of 0.1 to the corrected result. We show the corrected results in Fig. 15, where we include the results for aspect ratio N σ =N τ ¼ 6 directly at ðrTÞ min ¼ 0.5 as large open symbols. The screening mass is actually slightly above 2 times the perturbative Debye mass at NLO [Eq. (25)], indicated by the horizontal lines in the lower left corner for the scales μ ¼ πT, 2πT, and 4πT. We recall that in Fig. 15 we compare with 2 times the NLO Debye mass times a factor A ¼ 1. 25. This indicates that for F sub QQ we cannot clearly distinguish between the electrostatic and the asymptotic screening regimes numerically.
In order to test whether we have already reached the asymptotic behavior, we also include plane-averaged Polyakov loop correlators in the analysis, where the Polyakov loops are integrated over the xy plane.
Since the operators are separated only in the z direction, we have r ¼ z and rT ¼ zT in this case. Since the integration is only two-dimensional, the plane-averaged Polyakov loop correlator behaves as in the regime of asymptotic screening, where C 0 is a constant that diverges in the continuum limit. We did not subtract the asymptotic constant along the MD time history, as we found that it is not very helpful given numerical fluctuations in the large distance tail. Using the same approach as discussed for the point-to-point Polyakov loop correlation functions, we determine the screening mass from the plane-averaged Polyakov loop correlation function. Results for both operators are fairly consistent, although the plane-averaged correlator yields for small separations at low temperatures a systematically higher and at high temperatures a systematically lower screening mass by a few percent. These differences are of similar size as the errors in the smaller volume, as the difference between determinations in different volumes, and as well as our correction for the screening mass obtained at ðrTÞ min ¼ 0.25 [cf. the discussion preceding Eq. (D2)]. Therefore, we consider the screening masses as consistent and asymptotic for ðrTÞ min ≥ 0.5.
We show the local screening mass m S =T from F sub S in Fig. 33. Fits with Eq. (D1) as Ansatz yield m S < m D and χ 2 =d:o:f: > 1 for ðrTÞ min < 0.3. Both outcomes are consistent with the observation in Sec. IVA that screening is partially compensated by other medium effects at short distances. For larger ðrTÞ min , we always obtain χ 2 =d:o:f: < 1. m S =T for ðrTÞ min ∼ 0.3 − 0.6 is within errors fairly consistent with the perturbative Debye mass at NLO, as indicated by the three horizontal lines in the lower left corner of the figure representing m D for the scales μ ¼ πT, 2πT, and 4πT. Deviations can be attributed to the higher order terms discussed in Sec. IV B, which are missing in Eq. (D1). There is good numerical agreement between the screening masses for both volumes for ðrTÞ min ≲ 1. Even in the larger volume, we do not see an unambiguous plateau in m S =T for any of the temperatures, although m S =T is within errors consistent with a constant for ðrTÞ min ≳ 1. For this reason, we consider ðrTÞ min ≳ 1 as the asymptotic regime, although the onset of the asymptotic behavior seems to occur already somewhat earlier for lower temperatures. For ensembles with poorer signal-to-noise ratio or low temperatures (T < 2T c ), the fits may fail well before ðrTÞ min ∼ 1 and we cannot determine m S =T in the . We show the screening mass in units of the temperature, m S =T. The screening mass in the larger volume is shown with patterned bands, while the screening mass in the smaller volume is shown with filled symbols. The triplets of horizontal lines on the left side indicate the naive weak-coupling expectation m D =T for the scales μ ¼ πT, 2πT, and 4πT (solid, dotted, and dashed) using the NLO result [Eq. (25)]. The triplets of horizontal lines on the right side indicate the naive weak-coupling expectation for exchange of a bound state 2m D =T of two scalar modes, each with the mass m D =T. The scales are the same. Lastly, the horizontal solid bands indicate the estimate of the asymptotic screening mass obtained by shifting the screening mass for rT min ¼ 0.5 by a temperatureindependent constant. asymptotic regime. As for the case of the F sub QQ , we estimate a correction to the screening mass determined at shorter distances. We observe that the difference between m S =T determined for ðrTÞ min ¼ 0.5 and ðrTÞ min ¼ 1.3 is about 0.75. We use this difference to correct for screening masses m S =T determined at ðrTÞ min ¼ 0.5 on lattices with aspect ratio N σ =N τ ¼ 4 and assign a systematic uncertainty of 0.2 to the corrected result. The corrected results are shown as the solid bands in Fig. 33 and appear to underestimate the actual asymptotic screening mass somewhat at low temperatures, T ≲ 3T c . The screening masses in Fig. 14 are obtained with this procedure. We include the results for aspect ratio N σ =N τ ¼ 6 directly at ðrTÞ min ¼ 1.3 as large open symbols. We consider this estimate of the asymptotic screening mass as sufficiently accurate for T > 3T c .