Static $\bar{Q}Q$ pair free energy and screening masses from correlators of Polyakov loops: continuum extrapolated lattice results at the QCD physical point

We study the correlators of Polyakov loops, and the corresponding gauge invariant free energy of a static quark-antiquark pair in 2+1 flavor QCD at finite temperature. Our simulations were carried out on $N_t$ = 6, 8, 10, 12, 16 lattices using Symanzik improved gauge action and a stout improved staggered action with physical quark masses. The free energies calculated from the Polyakov loop correlators are extrapolated to the continuum limit. For the free energies we use a two step renormalization procedure that only uses data at finite temperature. We also measure correlators with definite Euclidean time reversal and charge conjugation symmetry to extract two different screening masses, one in the magnetic, and one in the electric sector, to distinguish two different correlation lengths in the full Polyakov loop correlator.


Introduction
At high temperatures strongly interacting matter undergoes a transition where colorless hadrons turn into a phase dominated by colored quarks and gluons, the quark gluon plasma (QGP). Recently, lattice simulations have shown that this transition is a crossover [1] and its characteristic temperature has also been determined [2][3][4][5][6]. Deconfinement properties of the transition can be studied by infinitely heavy, static test charges. At zero temperature a heavy quark and antiquark pair forms a bound state (quarkonium state), but above the deconfinement temperature, color screening and collisions with the medium would reduce the binding between the quark and the antiquark, eventually causing a dissociation. As proposed by Ref. [7], the behavior of quarkonia can signal deconfinement and QGP production in heavy ion experiments. Moreover, the different melting temperatures of the different states can be used as a thermometer, analogously to the spectral analysis of stellar media in astrophysics, where the absence and presence of the different spectral lines is used to determine the temperature.
In medium quarkonium properties are characterized by the corresponding spectral functions, studied in several works. However, extracting spectral functions from Euclidean meson correlators (i.e. the analytic continuation of the correlator to real time) is a difficult, ill-posed problem. Nevertheless, lattice studies of charmonium spectral functions using the Maximum Entropy Method have been carried out on numerous occasions [8][9][10][11][12][13][14][15][16][17][18][19]. A recent, detailed study of charmonium spectral functions in quenched QCD can be found in [13]. Results regarding spectral functions with 2 flavours of dynamical quarks can be found in Refs. [14,15]. A recent study of charmonium spectral functions in 2+1 flavour QCD is [16,20]. Bottomonium spectral functions have also been studied with the help of NRQCD [17][18][19].
Since the direct determination of the spectral function is difficult, one can study inmedium properties of quarkonium using approximate potential models. There are numerous proposals in the literature for lattice observables which can provide input to these models. The so-called singlet and octet potentials have been proposed [21][22][23][24][25], and studied on the lattice, but these are not gauge invariant, therefore extracting physical information from them is not straightforward. There was also a suggestion about using the analytic continuation of the Wilson-loop [26,27], that, however, has similar problems as the direct reconstruction of the spectral functions. Here, we calculate the gauge invariant static quark-antiquark pair free energy, a non-perturbatively well defined quantity, that carries information on the deconfinement properties of the QGP.
In the present paper we determine the free energy of a static quark-antiquark pair as a function of their distance at various temperatures. We accomplish it by measuring the Polyakov loop correlator [28], which gives the gauge invariantQQ free energy 1 as: (1.1) In the above formula, x runs over all the lattice spatial sites, and the Polyakov loop, L(x), is defined as the product of temporal link variables U 4 (x, x 4 ) ∈ SU (3) 2 : (1.2) or in the continuum formulation, as a path ordered exponential of the integral of the gauge fields: The leading order term to the correlator of Polyakov loops is a two gluon exchange diagram. It was first calculated at leading order in the dimensionally reduced effective theory (EQCD 3 ) [29]. Due to the two gluon exchange, the r dependence in leading order is exp(−2m D r)/r 2 , where m D is the Debye screening mass. This suggests that in the r → 0 limit, where perturbation theory is applicable, the correlator should behave as 1/r 2 . However, this is not the r → ∞ asymptotic behavior, which we need to fit the correlation length on the lattice. The reason is simple: even in the weak coupling limit, at distances larger than (g 2 T ) −1 the physics of magnetic screening becomes dominant. From the then applicable 3D effective pure Yang-Mills theory, Ref. [30] argued that at high temperature, the large distance behavior is exp(−m M r)/r, where m M is the magnetic screening mass. This was confirmed by 2 flavour lattice simulations (using a somewhat heavy pion) in [31]. 1 More precisely, the excess free energy that we get when inserting two static test charges in the medium. 2 In the literature, a factor of 1 Nc is often included in the definition. Including this factor leads to a term in the static quark free energy that is linear in temperature.
A related problem is that for the gluon self-energy, perturbation theory breaks down at the O(g 2 T ) order because of infrared divergences. This term contains contributions from magnetic gluons. Therefore, the perturbative definition of the screening mass, as a pole in the gluon propagator, is of limited use, since pertubation theory breaks down ( [32]). It is better to define the screening masses as inverse correlation lengths in appropriate Euclidean correlators. In order to investigate the effect of electric and magnetic gluons separately, one can use the symmetry of Euclidean time reflection [32], that we will call R. The crucial property of magnetic versus electric gluon fields A 4 and A i is that under this symmetry, one is intrinsically odd, while the other is even: Under this symmetry the Polyakov loop transforms as L R − → L † . One can easily define correlators that are even or odd under this symmetry, and thus receive contributions only from the magnetic or electric sector, respectively [31,32]: We can further decompose the Polyakov loop into C even and odd states, using A 4 and L C − → L * as: Next, we note that Tr L E+ = 0 = Tr L M − , so the decomposition of the Polyakov loop correlator to definite R and C symmetric operators contains two parts 3 . We define the magnetic correlation function as: 9) and the electric correlator as 4 : Then, from the exponential decay of these correlators, we can define the magnetic and electric screening masses. Note that with our definition TrL M + = Re TrL and TrL E− = i Im TrL , and: from which it trivially follows that if the magnetic mass screening mass is lower than the electric mass, we will have C(r, T ) − C(r → ∞, T ) asymptotic to C M + (r, T ) as r → ∞, or equivalently, the highest correlation length in C equal to that of C M + .
As for the asymptotic form of these correlators, similar arguments apply as with the full Polyakov loop correlator. In the high temperature limit the asymptotic behavior will be dominated by a glueball mass in the 3D effective Yang-Mills theory [30,32], but because of the symmetry properties, the quantum numbers carried by the glueballs will be different. We will therefore fit the ansatz: to extract screening masses, noting that the ansatz in principle is only motivated at high temperatures, where the effective field theory applies. Nevertheless we find that even close to T c the ansatz describes the large r tails of our lattice data well.

Simulation details
The simulations were performed by using the tree level Symanzik improved gauge, and stout-improved staggered fermion action, that was used in [33]. We worked with physical quark masses, and fixed them by reproducing the physical ratios m π /f K and m K /f K [33]. Compared to our previous investigations of Polyakov loop correlators, reported in the conference proceedings [34], here we used finer lattices, namely we carried out simulations on N t = 12 and 16 lattices as well as on N t = 6, 8, 10 lattices. Our results were obtained in the temperature range 150 MeV ≤ T ≤ 450 MeV. We use the same configurations as in Ref. [5] and [35]. Figure 1 summarizes our statistics.

Renormalization procedure and continuum extrapolation
After measuring the Polyakov loop correlator C(r, T ) at T = 1/(N t a) temperature, we computed the unrenormalized free energy according to FQ Q = −T ln C(r, T ). The a(β) function was taken from the line of constant physics, along which we kept the ratios of the physical values of m π , f K and m K fixed at zero temperature. Detailed description of the determination of the line of constant physics can be found in Ref. [35].
Approaching the continuum limit, the value of the unrenormalized free energy diverges. In order to eliminate the divergent part of the free energy renormalization is needed. There are various proposals in the literature for this renormalization procedure. Earlier works [22][23][24] matched the short distance behavior to the T = 0 static potential, but this is ambiguous. A more precise definition is to require that the T = 0 potential vanishes at some distance [3,34]. This would require a precise determination of the potential at T = 0. Here, though, we use a renormalization procedure based entirely on our T > 0 data, similarly to Ref. [36]. The data contains a temperature independent divergent part from the ground state energy. The difference between the value of free energies at different temperatures is free of divergences. Accordingly, we define the renormalized free energy as: with a fixed T 0 . This renormalization prescription corresponds to the choice that the free energy at large distances goes to zero at T 0 , and is implemented in two steps. In the first step we have: where the one quark free energy F Q satisfies: Note, that this first step of the renormalization procedure is completely straightfoward to implement, at each simulation point in N t and β we just subtract the asymptotic value of the correlator. This gives us a UV finite quantityFQ Q , however we don't call this the renormalized free energy, since compared to equation (3.1) it contains less information. Namely at all temperaturesFQ Q (r = ∞, T ) = 0. The correlation length is the same as with definition (3.1), but the information of the asymptotic value (that is the single heavy quark free energy) is lost. That information however is retained in the second step: where the renormalized one heavy quark free energy is:  Doing the renormalization in two steps like this has a technical reason that will be explained shortly.
Let us first mention that this Polyakov loop correlator behaves similarly to the baryon correlators in imaginary time do: at large values of r we can get negative values of C at some configurations 5 . For this reason, it is highly desirable to use gauge field smearing which makes for a much better behavior at large r, at the expense of unphysical behavior at small r. For this reason, we measured the correlators both without and with HYP smearing. We expect that outside the smearing range (i.e. r ≥ 2a) the two correlators coincide. This is supported by Figure 2. Therefore we use the smeared correlators for r ≥ 2a and the unsmeared ones for r < 2a.

Single heavy quark free energy
First, we discuss the implementation of the renormalization of the single heavy quark free energy, equation (3.5). Notice that if we implemented the renormalization condition (3.1) directly, then we would just need to subtract 2F Q (β, T 0 ) from the unrenormalized free energy, so at first sight it looks like we are doing some unnecessary rounds by doing this in two steps. What we gain by this is that we can extend the temperature range, at which we can do the continuum limit. To understand this statement let us look at  corresponding to the green line in the figure, this gives us 5 points from the curve F Q (β, T 0 ). According to equation (3.5) this is what we have to subtract from the bare free energy at this value of β to get the renormalized single quark free energy. The disadvantage of the green curve, is that the β range it covers is rather limited. So, if we want to be able to make a continuum limit from say the N t = 8, 10, 12 lattices, the temperature range we can cover is rather limited as well. The lowest temperature we will be able to do a continuum limit at will be (6/8) × 200MeV = 150MeV, and the highest temperature will be (16/12) × 200MeV = 266MeV.To do a continuum limit at higher temperatures, we need the F Q (β, T 0 ) curve at higher values of β, and at first, it looks as like that would need runs at higher values of N t . This is not feasible, but there is a simple tricj to extend the temperature range. Clearly, if we call the continuum limit of the single quark free energy than, for any value of T : is just a number 6 . We can use this fact to extend the temperature range of the continuum limit by using different values of T 0 , that is different renormalization prescriptions, and shift them together by the value of the RHS of equ. (3.7). This is the procedure that we will follow.
To implement equation (3.5), we first calculate F Q (β, N t ) or equivalently FQ Q (r → ∞, β, N t ) from equation (3.3). Then at each N t we interpolate to the β value corresponding to the temperature T 0 , giving us some points of the function F Q (β, T 0 ). Finally, we interpolate these F Q (β, T 0 ) points in β, obtaining the final curve we can use for the renormalization. This procedure is illustrated on Figure 3 (left). When doing this interpolation we take into account the error on the data points of F Q (β, N t ) via the jacknife method. The statistical errors of the single quark free energy are very small, meaning that the interpolation method gives a comparable error to the final interpolated value. We estimate the systematic error of the interpolations by constructing different interpolations. For interpolations of the F Q (β, N t ) curves we use linear and cubic spline interpolations (for each value of N t ), and for the interpolation of F Q (β, T 0 ) we use different polynomial interpolations(order 1,2), cubic spline and barycentric rational function interpolation. In total this means 2 5 × 4 = 128 different interpolations, than for interpolating the bare F Q we use spline and linear interpolations, so for the final renormalized values we have in total 128×2 = 256 different interpolations. All interpolations are taken to have the same weight. We use the median of this as the estimate, and the symmetric median centered 68% as the 1 systematic error estimate [37]. The statistical and systematic errors turn out to be of the same order, and are than added in quadrature.
After doing this procedure, the β range in which we can interpolate the F Q (β, T 0 ) curve is limited, therefore, the temperature range where we can do the continuum extrapolation is limited. To extend the temperature range where we can calculate the single heavy quark free energy, we use the fact that the single heavy quark free energies at different temperatures differ only by an additive constant in the continuum. Therefore we use different values of T 0 to do the continuum extrapolation, and shift all those curves to the position of the 200MeV curve. We used 5 different values of T 0 , namely, 170MeV, 200MeV, 240MeV, 320MeV, and 390MeV. The results of this analysis can be found in Figure 3 (right).
For the continuum limits, we use the N t = 8, 10, 12 lattices, that are available at all temperatures. We use the N t = 16 lattice to estimate the systematic error of the continuum extrapolation, where it is available. If: Finally, we mention that the determination of the continuum limit of the Polyakov loop, or equivalently, that single static quark free energy is already available in the literature. For two recent determinations of the Polyakov loop see Refs. [5,38]. The difference is that here we take the continuum limit at significantly higher temperatures.

HeavyQQ pair free energy
Next, we turn to the determination ofFQ Q defined in equation 3.2. This quantity is UV finite, and goes to 0 as r → ∞. Similarly to the single quark free energy, the determination ofFQ Q at a given value of T and r requires two interpolations. At first we are givenFQ Q at several values of T , at each T we have a different value of the lattice spacing. If we want to know the value ofFQ Q at (T, r) = (T * , r * ) at some value of N t , first we do an interpolation in the r direction to the value r * at each given T , then we do an interpolation in the T direction, where the node points for the interpolations are the interpolants in the previous step. The statistical error than can be estimated by constructing these interpolations to every jacknife sample. For systematic error estimation we try different interpolations in the r and T directions. In the r direction we have: polynomials of order 1,2,3,...,7 and a cubic spline, in the T direction we have polynomials of order 1,2,3 and cubic spline. This is in total 4 × 8 = 32 different interpolations. Just as before, we use the median of these values as the estimate, and the symmetric median centered 68% as the 1σ systematic error estimate. Like in the case of the single heavy quark free energies, the statistical and systematic errors turn out to be of the same order, and are then added in quadrature.
Next, we do the continuum extrapolation. Here we also take a similar approach as in the previous subsection. For the continuum extrapolations, we use the N t = 8, 10, 12 lattices, that are available at all temperatures. We use the N t = 16 lattice to estimate the systematic error of the continuum extrapolation, exactly like before. Also, where the N t = 16 lattices are not available, we estimate the systematic error, as in the previous section, by the average of the systematic error at the points where we do have N t = 16 lattices (approximately 7%). The linear fits of the continuum limit extrapolations all have good values of χ 2 .
Next, we add the values of 2F Q , determined in the previous subsection, and visible in Figure 3 to the free energy values to obtain the final results in Figure 4 (errors are added in quadrature). Note, that the N t = 6 lattices were only used in the whole analysis to extend the β range of the renormalization condition for the single quark free energy.

Magnetic and electric screening masses
We continue with the discussion of the electric and magnetic screening masses obtained from the correlators (1.9) and (1.10). For this analysis we only use lattices above the (pseudo)critical temperature, since that is the physically interesting range for screening. Next, we mention that for this analysis, we only use the data with HYP smearing, since we are especially interested in the large r behavior. Before going on to the actual fitting  . Continuum values of the staticQQ free energy at different temperatures. Note that the curves seem to tend to the same curve as r → 0, corresponding to the expectation that UV physics is temperature independent.
procedure of the screening masses let us first illustrate some simple relations, with the raw lattice data of the electric and magnetic correlators. First C E− (r, T ) C M + (r, T ) as r → ∞, or equivalently, that the electric screening mass is larger than the magnetic one. This can be seen on Figure 5. The next observation is that the screening masses in both channels are approximately proportional to the temperature. This can be seen on Figure  6. Both of these facts are expected to hold at high temperatures, but these lattice results suggest that they hold at lower temperatures as well.
Next, we turn to the actual determination of the screening masses. So far there has been one determination of electric and magnetic screening masses on the lattice using the non-perturbative definition given by ref. [32]. That study used 2 flavours of Wilson fermions with a somewhat heavy pion, and did not attempt a continuum extrapolation [31].
Since the masses are expected to be proportional to the temperature, the natural distance unit in this problem is rT , so we give limits on the range of the fits in these units. For the correct determination of the screening masses, special care is needed in the choice of the fit interval. To find the proper minimum rT value of the fits, we use hypothesis testing, similar to that in Ref. [39]. If the fits are good, than the value of χ 2 , defined as: should have a χ 2 distribution, with the appropriate degrees of freedom. Here C ij is the Since the x axis of this plot is rT , if one assumes a Yukawa form of the correlator, than the slopes of these curves are just m M /T and m E /T respectively. The fact that the graphs are approximately parallel straight lines suggests that these ansatzes are approximately correct, and that the masses are approximately proportional to the temperature.
covariance matrix. In this case the quantity Probability density of χ 2 (x)dx, (4.2) should have a uniform distribution on [0, 1]. If we fix the range of all the fits in rT units, each fit (at some value of N t and β) gives one pick from a supposed uniform distribution in Q. This is equivalent to having multiple picks from the same uniform distribution. We will test this hypothesis with a Kolmogorov-Smirnov test for the uniform distribution. Here one determines the maximum value of the absolute difference between the expected and measured cumulative probability distributions. This is then used to define a significance level or probability that the measured distribution can indeed be one originating from the expected uniform distribution. These probabilities are listed in Table 1. We will only use value of (rT ) min where the Kolmogorov probability is at least 0.3. This test tells us, that for systematic error estimation, we will have, for the magnetic correlator (rT ) min going from 0.465 to 0.61, and for the electric correlator we have (rT ) min going from 0.35 to 0.43 7 .
At this point we mention that for the continuum limit we will not use the N t = 16 lattices, because the mass fits there have huge error bars. Nevertheless, when the continuum limit is done, we will see that the values of the masses at the N t = 16 lattices are consistent with the continuum estimates. Also, if we use them, we get the same results, because they do not give a contribution to the continuum limit, due to the big errors. Now that we have estimated the proper rT range of the fits, we go on to the fitting of the masses. The results of the fits at different values of N t can be seen in Figure 7. The systematic errors come from changing the lower limit of the fit, in the case of the magnetic correlator, from (rT ) min = 0.465 to (rT ) min = 0.61, and in the case of the electric correlator, from (rT ) min = 0.35 to (rT ) min = 0.43. The results coming from different values of (rT ) min are weighted using the Akaike Information Criterion(AIC) [40]. The median of the weighted histogram gives the central value, and the central 68% the systematic error estimate. Note that using the Q values as weights or uniform weights gives a very similar 7 (rT )max was fixed on both cases. Increasing (rT )max results in a less precise covariance matrix and correspondingly, somewhat worse χ 2 values, but consistent screening masses. For example, if for the magnetic correlator we choose (rT )max = 1 instead of 0.9, the final value of the Kolmogorov-Smirnov probability in Table 1 will not be 96%, but 38% instead. Nevertheless the growing trend in the probabilities will be the same. Also, we will get the same results within uncertainties. result. The statistical error comes from a jacknife analysis with 20 jacknife samples. The two errors turn out to be of similar magnitude (with the statistical error being somewhat bigger) and are then added in quadrature.
Next, we fit linear functions to all screening masses at all values of N t , and use these to do a continuum extrapolation from the N t = 8, 10, 12 lattices. Taking into account the errors of the linear fits, all χ 2 values of the continuum limits are very good. The continuum limit, in addition to the statistical error, also has a systematic error estimated, from doing a 2 point linear extrapolation from the N t = 12, 10 lattices, and taking the difference of the extrapolated value from fitted value to the N t = 8, 10, 12 lattices 8 . The statistical and systematic errors are added in quadrature.
We finish this section by comparing our results to those from earlier approximations in the literature. For comparison let us use our results at T = 300MeV ≈ 2T c . Here we have: This work 3D EFT, N f =3 AdS/CFT WHOT QCD, N f =2, N t =4, m π /m ρ =0.65 Figure 8. The continuum extrapolations of the screening masses and the ratio of the screening masses. For the ratio m E /m M we also included different estimates from the literature: Lattice results from Ref. [31], dimensionally reduced 3D effective field theory results from Ref. [41], and results from N = 4 SYM plasma with AdS/CFT from Ref. [42].

Conclusions
In this paper we have determined the renormalized static quark-antiquark free energies in the continuum limit. We introduced a two step renormalization procedure using only the finite temperature results. The low radius part of the free energies tended to the same curve, corresponding to the expectation that at small distances, the physics is temperature independent. We also calculated the magnetic and electric screening masses, from the real and imaginary parts of the Polyakov loop respectively. As expected, both of these masses approximately scale with the temperature as m ∝ T , with m M < m E , therefore, magnetic contributions dominating at high distances. The values we got for the screening masses are close to the values from dimensionally reduced effective field theory.