Un-screened forces in Quark-Gluon Plasma ?

,


I. INTRODUCTION
At very high temperatures the strongly interacting matter undergoes a transition to a new state called quark-gluon plasma (QGP).Creating and studying the properties of QGP is the goal of large experimental programs in heavy-ion collisions at RHIC and LHC [1].
The question of in-medium modifications of the forces between heavy quark Q and antiquark Q generated a lot of interest since the seminal paper by Matsui and Satz [2].They conjectured that color screening in QGP will make the Q Q interaction short ranged, and therefore quarkonium states cannot be formed in QGP.Thus, QGP formation in heavy-ion collision will lead to quarkonium suppression.The study of quarkonium production in heavy-ion collisions is a large part of the experimental heavy-ion program, see e.g.Ref. [3] for a recent review.
The idea of having a screened potential between heavy quarks in QGP is closely related to the exponential screening of the free energy of infinitely heavy quarks in QGP, which is well established by lattice QCD calculations, see e.g.Ref. [4] for a review.However, the free energy of heavy quarks describes the in-medium interaction of heavy quarks at macroscopic time scales much larger than the inverse temperature.For understanding the quarkonium properties in QGP one needs to know if and how the heavy Q Q potential is modified at scales comparable to the internal time scale of quarkonium.The effective field theory (EFT) approach provides a natural framework to address this problem at high temperatures when the weak-coupling approach is applicable [5,6].Depending on the separation of the bound-state scales and the thermal scales the heavy Q Q potential can be modified by QGP and also acquire an imaginary part.In general, however, the real part of this potential does not have a screened form in this approach [6].How to study the modification of heavy Q Q interactions in QGP beyond weak coupling remains an unsolved problem.However, we could define the heavy Q Q potential at non-zero temperature (T > 0) in analogy with the zero temperature (T = 0) case in terms of the Wilson loops of size τ × r [7].We can write the following spectral representation of the Wilson loops in terms of the r-dependent spectral function (1) The distance r between the heavy Q and Q acts as the label of the spectral function.At T = 0, the spectral function's lowest delta function peak corresponds to the ground state potential.We expect that there is a dominant, broadened peak in the spectral function for not too high temperatures; its position and width determine the real and imaginary parts of the potential, respectively [7].For very high temperatures the spectral function may lack a well-defined peak such that a potential cannot be defined.While the relation between the above defined complex potential and the EFT concept of the complex potential is an unsolved problem, too, the existence of a well-defined peak in ρ r (ω, T ) is necessary, yet not a sufficient condition for a potential picture of heavy quarkonium at T > 0.
In this paper we present calculations of the real part of the potential at T > 0 in 2+1 flavor QCD using the lattice QCD approach and estimate the imaginary part.There have been several attempts to calculate the complex potential at T > 0 both in quenched QCD [7,8] as well as in 2+1 flavor QCD [9,10].The state of the art calculation of the complex potential in 2+1 flavor QCD has been performed using lattices with temporal extent N τ = 12, and thus at a single lattice spacing per temperature.The new results are based on several lattice spacings and several values of N τ in the range N τ = 16 − 36.The rest of the paper is orgnized as follows.In section II we give some details of the lattice QCD calculations.Section III presents the analysis of the lattice results and the main results of the study, while section IV contains our conclusions.Many technical details of the calculations are discussed in the Appendix.

II. DETAILS OF THE LATTICE QCD CALCULATIONS
In lattice QCD one often considers correlators of Wilson lines in Coulomb gauge instead of Wilson loops since these contain the same physical information and are less noisy, see the discussions in Ref. [10] and Appendix A. We calculated the Wilson line correlators W (τ, r, T ) in 2+1 flavor QCD using highly improved staggered quark (HISQ) action [11] and Lüscher-Weisz action [12,13]  MeV in the continuum limit.Furthermore, the calculations have been performed for three different lattice spacings corresponding to the following values of bare lattice gauge coupling β = 10/g 2 0 = 7.596, 7.825 and 8.249.The lattice spacing and thus the temperature scale T = 1 /(aNτ ) has been fixed using the r 1 -scale determined in Ref. [14] with the value r 1 = 0.3106 fm obtained in Ref. [15].The value of the strange quark mass was obtained from the parametrization of the line of constant physics from Ref. [16].We this scale setting for the lattice spacing we obtain: a(β = 8.249) = 0.0280 fm, a(β = 7.285) = 0.0404 fm and a(β = 7.596) = 0.0493 fm.For the finest lattices the spatial size of the lattice is N s = 96, while for the two coarser lattices we use N s = 64.The temporal size of the lattice is varied in the range N τ = 16 − 36, which corresponds to the temperature range 153 MeV ≤ T ≤ 352 MeV.Further details about the parameters of the lattice calculations are given in Appendix A.
For the smallest lattice spacing, we only consider T ≥ 195 MeV, i.e. temperatures well above the chiral crossover temperatures.For these temperatures we do not expect significant quark mass dependence of the Wilson line correlators.Therefore, the calculations for the smallest lattice spacing have been performed only with m l = m s /5, while for the coarser lattices we use m l = m s /20.As discussed later, we do not see any m l dependence of the Wilson line correlator for T ≥ 195 MeV, as expected.
For lattices with large temporal extents, employed in this study, noise reduction methods have to be used.We use gradient flow [17] for noise reduction.To reduce noise even further, we require that W (τ, r, T ) at large r/a is a smooth function of r/a, and replace it for each value τ by a corresponding local r/a interpolation.We verified that this procedure does not introduce bias in our analysis by varying the interpolation intervals and comparing to the results that do not use interpolation.Further details on the noise reduction techniques are presented Appendix A.
To aid the reconstruction of the spectral function we also performed calculations on N τ = 64 and N τ = 56 lattices, which we refer to as T = 0 lattices.

III. ANALYSIS AND RESULTS
To analyze the Wilson line correlator W (τ, r, T ) in Eq. (1) it is useful to consider the effective mass defined as where the last equation applies to the case of nonzero lattice spacing.At T = 0, the effective mass decreases with increasing τ , and reaches a plateau for sufficiently large τ , since the spectral function is positive definite and has the lowest ground state delta function peak followed by many excited states for ω above the ground state.We show the results for the effective masses in Fig. 1.We see that at T = 0 the effective mass decreases with increasing τ , with the exception of the data at smallest τ , and approaches a plateau for τ around 0.5 fm.The non-monotonic behavior is due to the smearing artifacts coming from the gradient flow, as discussed in Appendix B. Except for very small τ , m eff decreases at T > 0 with increasing τ for all τ values and does not reach a plateau.This means that there is no stable ground state at non-zero temperature.We see from Fig. 1 that the effective masses show neither lattice spacing nor sea quark mass dependence for T > 200 MeV.This implies that for these temperatures using m l = m s /5 is equivalent to using the physical light quark mass and that our results are essentially in the continuum limit.We also compared the effective masses corresponding to different lattice spacings at lower temperatures and found no dependence on the lattice spacing.At small τ the difference between the T = 0 or T > 0 effective masses is the smallest, and their τdependence is rather similar, see Fig. 1.Thus, we aim to constrain the corresponding contributions at T > 0 by using the T = 0 results.
Our objective is to extract information on a dominant peak in the spectral function corresponding to W (τ, r, T ) at T > 0. We choose an Ansatz [10] for the spectral function as where ρ high r (ω) is assumed to be a temperatureindependent part dominating at large ω.ρ peak r (ω, T ) describes a dominant peak encoding the complex potential at T > 0, while ρ low r (ω, T ) is a small, medium-dependent contribution below the dominant peak.
Fixing ρ high r (ω) to its T = 0 value effectively means subtracting it from the T > 0 result.We define the subtracted correlator as follows, which is solely determined by the mediumdependent contributions to the spectral function, (5) Since the T = 0 spectral function has the form we define via ) Thus, it is straightforward to estimate W high (τ, r) using single-exponential fits for A r and V (r, T = 0).The task of constraining ρ peak r (ω, T ) and ρ low r (ω, T ) is now reduced to the analysis of the τ -dependence of W sub (τ, r, T ).
The effective masses from W sub (τ, r, T ) at T > 0 are also shown in Fig. 1.The uncertainties in the effective masses due to the errors in the determination of the ground state contribution at T = 0 have been taken into account by combining these uncertainties with the statistical errors of the T > 0 calculations.We note that the non-monotonic behavior at small τ due to smearing artifacts is absent in these subtracted effective masses m sub eff (τ, r, T ), and therefore, these artifacts do not affect ρ peak r (ω, T ), see Appendix C. As discussed in Appendix C, m sub eff (τ, r, T ) would decrease linearly in τ if the ground state peak had a Gaussian shape [10].m sub eff (τ, r, T ) shows linear behavior in τ at small τ , indicating that the dominant ground state peak has broadened.Here we note, that the behavior of the subtracted effective masses obtained from the Wilson line correlators and from the Wilson loops is the same [10].
As discussed in Ref. [10] ρ low r (ω, T ) is the contribution to the spectral function at T > 0 which has support for energies well below the dominant peak and representing a heavy Q Q state propagating forward in Euclidean time interacting with a backward propagating light state from the medium.This contribution is much smaller than ρ peak r (ω, T ) but dominates the correlator at τ around 1/T .This part of the spectral function explains the rapid drop of m eff (τ, r, T ) at large τ [10] that can be seen in Fig. 1.
A physically appealing parametrization of ρ peak (ω, T ) is a Lorentzian form.However, a Lorentzian form is only valid in the vicinity of the peak.In general, we can assume that the correlator has a pole at some complex ω, so For ω ≃ ReV (r, T ) we can approximate Γ(ω, r, T ) by a constant: Γ(ω, r, T ) ≃ Γ L (r, T ).However, for ω values far away from the peak Γ(ω, r, T ) must quickly go to zero.The selfconsistent T -matrix calculation of heavy Q Q propagators indeed shows an exponential decrease of Γ(ω, r, T ) away from the peak [18].To incorporate this feature of the spectral function in our analysis we assume that ρ peak r (ω, T ) is given by Γ ) and is zero otherwise.Such a cut Lorentzian form gives rise to an almost linear behavior of m sub eff (τ, r, T ) at small τ , too, as required by the lattice data.The most general parametrization of ρ low r (ω, T ) would be a sum of delta functions at ω well below the dominant peak position.However, to describe our effective mass data even a single delta function at sufficiently small ω, ρ low r (ω, T ) = c low r (T )δ(ω − ω low r (T )) turns out as sufficient.With these forms of ρ peak r (ω, T ) and ρ low r (ω, T ) we fitted the lattice data on m sub eff (τ, r, T ) and determined the fit parameters ReV (r, T ), Γ L (r, T ), c low r (T )/A r (T ) and ω low r (T ).A sample fit is shown in Fig. 1 and details of the fits are discussed in Appendix C. We typically find that c low r (T )/A r (T ) < 5 • 10 −4 and decreases with decreasing r, while ω low r (T ) is between (1.8-3.8)GeV below the peak position ω = ReV (r, T ).
The results for ReV (r, T ) are shown in Fig. 2 indicating a temperature-independent real part in good agreement with the T = 0 potential.This is not completely unexpected, as m eff (τ, r, T ) at small τ is close to the vacuum result, c.f. Fig. 1.The peak position is insensitive to the detailed shape of ρ peak r (ω, T ), i.e. for a Gaussian form we find the same peak position within errors.Thus our lattice QCD results show that the potential's real part is unscreened.This observation supersedes conclusions drawn earlier by applying the Bayesian Reconstruction (BR) method [19] to older lattice data [20] with much larger statistical errors and larger discretization artifacts.There are distortions in ReV (r, T ) at the two shortest distances in lattice units (r = a, 2a), but these distortions are the same both at T = 0 or T > 0, see the discussion in Appendix B and Appendix C, and do no affect our conclusion about the ab-sence of screening.
As discussed above the imaginary part of the potential is defined as the width of the ground state peak at T > 0. If we knew the spectral function exactly we could fit it in the peak's vicinity with a Lorentzian form, whose width parameter would give the potential's imaginary part.This has been explicitly checked for the spectral function of an infinitely heavy Q Q pair calculated in hard thermal loop perturbation theory [21].Yet the correlator is sensitive to all parts of the spectral function, in particular to ρ low r (ω, T ) and to the tails of ρ peak r (ω, T ).For this reason, the parameter Γ L cannot be considered as ImV (r, T ).A better way to characterize ImV (r, T ) is to consider the cumulants of ρ peak r (ω, T ).The first two cumulants are defined as c 1 = ⟨ω⟩ and c 2 = ⟨ω 2 ⟩ − ⟨ω⟩ 2 , where ⟨. . .⟩ = dω . . . .In the case of the Gaussian, the second cumulant of the spectral function is the square of the width parameter.In the case of the cut Lorentzian, it is proportional to the square of the parameter Γ L .Furthermore, if c low r /A r is very small, ρ peak r (ω, T ) determines the behavior of the Wilson line correlator around τ = 0. Therefore, the second cumulant of ρ peak r (ω, T ) determines the slope of m sub eff (τ, r, T ) at small τ , which is well defined from the lattice data, see the Appendix D. Thus the square root of the second cumulant of ρ peak r (ω, T ) is a good proxy for the r and temperature dependence of ImV (r, T ).In Fig. 3 we show this proxy for ImV (r, T ) as a function of distance r for different temperatures.We scaled the x-and y-axes by the temperature in the two middle and right panels of Fig. 3.We see that for 180 MeV < T ≤ 352 MeV the numerical results for ImV (r, T ) scale with the temperature, i.e. the imaginary part of the potential depends only on rT and is proportional to the temperature.This is in qualitative agreement with the weak-coupling results.Since for rT ≃ 1 the imaginary part of the potential is larger than the temperature, the forces between heavy quarks are damped very quickly, i.e. on the time scale comparable to or shorter than the thermal scale.During that short time scale, the chromo-electric field between the heavy Q and Q cannot adjust itself to the medium.The chromoelectric force between the heavy quarks is simply damped away, and the heavy Q and Q will not interact.This picture of quarkonium melting is very different from the one proposed by Matsui and Satz.While ImV is quite large we still think the Q Q energy is well defined in the considered temperature interval.For if there would be no well defined dominant peak in the spectral functions, different static Q Q correlators would have quite different τ dependence.However, as shown in our previous study [10] this is not the case.

IV. CONCLUSION
We studied the complex heavy quark-antiquark potential at non-zero temperature in 2+1 flavor QCD using lattice calculations with a large temporal extent.We have found that contrary to some common expectations the real part of the potential is not screened for temperatures 153 MeV ≤ T ≤ 352 MeV.We also found that the dissipative effects on the chromo-electric forces between the heavy quarks, encoded in the imaginary part of the potential are very large and likely will lead to quarkonium dissolution.
As already mentioned in the introduction, the lack of screening in the real part of the potential is expected in weak-coupling limit for rT < 1 [6].Our study shows that this also holds nonperturbatively.Furthermore, a numerical evaluation of the weak-coupling result of the thermal correction to the real part of the potential shows that this correction is quite small.For T > 500 MeV and rT < 0.4, where the weak coupling result of Ref. [6] may be applicable, the thermal correction to the real part of the potential is smaller than 0.1T .In this appendix we discuss further details of our lattice QCD calculations.The parameters of the lattice calculations including the lattice volume and the quark masses are given in Tables I, II, and III.The gauge configurations used in this study have been generated using a rational hybrid Monte-Carlo algorithm [23] with grants from PRACE on Juwels Booster and Marconi 100 and NERSC on Perlmutter using the SIMULATeQCD code [22].We also used the MILC code on Cori at NERSC to generate the gauge configurations.Some of the gauge configurations have been generated on the USQCD cluster in JLab.After removing the initial trajectories for thermalization we arrived at the data set in Tables I, II, and III.Every 5th trajectory has been used for N σ = 96 and every 10th trajectory for N σ = 64.
On the generated gauge configurations we calculated Wilson line correlators in Coulomb gauge with the aim of determining the static quarkantiquark (Q Q) potential.We use Wilson line correlators instead of Wilson loops because these are much less noisy and provide more convenient access to distances at non-integer multiples of the lattice spacing.At T = 0 both Wilson loops and Wilson line correlators in Coulomb gauge have been used for the determination of the Q Q potential, see e.g.Refs.[14,[24][25][26][27].In the case of Wilson loops, smearing should be applied to the spatial gauge links entering the Wilson loops in order to obtain a reasonable signal.In Ref. [10] both Wilson lines and Wilson loops with threedimensional hyper-cubic (HYP) smearing [28] in the spatial gauge links have been studied at nonzero temperature.It was found there that the behavior of the Wilson line correlators and Wilson loops is fairly similar except for small τ , where sensitivity to excited states is different, similar to the T = 0 case [29].At T > 0 there are also some differences between the behavior of Wilson loops and Wilson line correlators at τ ≃ 1/T , which are, however, not related to Q Q potential as discussed below.Thus both Wilson lines in Coulomb gauge and Wilson loops encode the same temperature modification of the Q Q potential.In Ref. [10] the calculations of the Wilson lines have been performed on N τ = 12 lattices.
Since we use much larger N τ in this study, also the temporal links have to be smeared.We use gradient flow [17] for the smearing of the temporal gauge links.More precisely we use Zeuthen flow [30].For flow time τ F the gauge links are smeared in a radius √ 8τ F .This radius should be much smaller than the inverse temperature.We use different flow times corresponding to the flow radius in the range a − 2.53a and study the sensitivity of our results to the flow time.For the final results presented in the paper, we use the smallest flow time that gives an acceptable signal.Since the signal deteriorates with increasing N τ we use larger flow time for large N τ .The range of flow times and the specific values of flow times for which we show the final result are presented in Tables I, II, and III for β = 8.249, 7.825 and 7.596, respectively.
After performing the gradient flow we fix the Coulomb gauge.The precision of Coulomb gauge fixing was set to 10 −6 .We also note that neither is the gradient flow the only option to smear the temporal gauge links nor is it a problem to fix the Coulomb gauge before performing the gradient flow, when studying Wilson line correlators.Previously we used HYP smearing after gauge fixing for the temporal gauge links when calculating the Wilson line correlators at T > 0 [31] and found that the temperature and the τ dependence of the correlators are similar to that reported here.Thus even though smearing destroys the gauge fixing condition to some extent, the qualitative behavior of the Wilson line correlators is not affected.This implies that our findings are neither sensitive to the details of gauge link smearing nor to details of the Coulomb gauge fixing.useful to consider the effective masses defined in Eq. (2.The Wilson line correlators require multiplicative renormalization which corresponds to an additive normalization of the effective masses that is proportional to 1/a.This normalization can be fixed by requiring for each lattice spacing that the Q Q potential at T = 0 is equal to a prescribed value for one given distance.Here we use the prescription V (r = r 0 ) = 0.954/r 0 , where r 0 is the Sommer scale, which for 2+1 flavor QCD is r 0 = 0.468(4) fm [26].This normalization condition was used in our previous studies [16,26,32].The normalization constant depends on the amount of smearing, i.e. the coefficient 2c Q of the 1/a divergence is smearing dependent.The larger the amount of smearing, the smaller the coefficient of the 1/a divergence becomes.For unsmeared Wilson line correlators the coefficient c Q was determined in Ref. [32] for several beta values including, the two lowest ones used here, namely c Q (β = 7.596) = 0.3545 (11) and c Q (β = 7.825) = 0.3403 (12).Interpolating the results for c Q from Ref. [32] with cubic polynomial we estimate c Q (β = 8.249) = 0.3144 (10).In Fig. 4 (top) we show the un-renormalized effective masses at T = 0 for β = 8.249 at different flow times.The improvement in the signal with increasing flow time at large τ is obvious from the figure.We also see that the effective masses decrease with increasing flow time as one would expect based on the discussions above.There is a non-monotonic behavior of the effective masses in τ for τ /a = 1 − 3.This is due to the fact that the gradient flow distorts short distance physics and potentially can lead to non-positive definite spectral function for very large ω.However, for not too large ω there is no sign of positivity violation in the spectral function since the effective masses approach plateaus from above for τ /a > 3.This means that the gradient flow does not lead to artifacts in the determination of the Q Q potential at T = 0.With a constant, flow-time dependent shift the effective masses for different τ F can be collapsed to one line, except for very small τ , where there are τ F -dependent distortions due to  gradient flow.This is demonstrated in Fig. 4 (bottom).We determine this shift by fitting the difference in the effective masses calculated at different flow times to a constant for τ /a = 7−18 for β = 8.249 and τ /a = 7 − 15 for the two smaller values of β.This constant shift should amount to the difference in the additive normalization of the Q Q potential, and therefore, should be independent of Q Q separation, r, apart from the distortions at small r due to smearing.In Fig. 5 we show the relative shifts as a function of r for β = 8.249 .We see that for very small r there is some dependence on the value of r implying that there are distortions in the zero temperature potential at these distances due to smearing as expected.Namely, when τ F /a 2 ≤ 0.2 we see distortion for r/a < 2, while for larger flow time we see distortions for r/a < 3. To demonstrate this in Fig. 6 we show the zero temperature potential for different smearing levels for relatively small r values.We see from the figure that that except for the smallest distance the potential does not depend on the smearing level including the case of no smearing.We found that the situation for the other two β values is the same.
In addition to the gradient flow, we use polyno- mial interpolations to reduce fluctuations in the Wilson line correlators.For fixed τ the Wilson line correlators should be a smooth function of r apart from the effects of breaking of rotational symmetry on the lattice.For Symanzik gauge action these effects are smaller than the statistical errors for r/a > 3 [29,33].Therefore, it is natural to require that the data on the Wilson line correlators are smooth functions of r at a fixed value of τ .By imposing this requirement we effectively reduce the fluctuations in the original data set since nearby r values usually correspond to very different path geometries and thereby suffer from quite independent gauge noise.We perform second order polynomial interpolations in a limited range of distances, ∆r around a target value of r and replace the original datum with the interpolated value.We take into account that, with increasing distance, there are many different separations that are close to the target value of r and adjust ∆r as we vary r.This additional noise reduction and the interpolation procedure are demonstrated in Fig. 7.In fact the result on the effective masses shown in Fig. 4 also incorporate the noise reduction from the interpolations.Because of the use of the above noise reduction the determination of the Q Q potential at zero temperature is now more accurate.Therefore, we re-calibrated the central value of the constant c Q and used the following values in the present analysis: c Q (β = 7.596) = 0.3552, c Q (β = 7.825) = 0.3401 and c Q (β = 8.249) = 0.3135.These values agree with the one quoted above within errors.
To check that interpolations do not introduced additional bias we performed the analysis by doubling the interpolation range in r, and also obtained the zero temperature potential without any interpolations.The results are shown in Fig. 8.As one can see the zero temperature potential is not sensitive to these changes.Doubling the interval in the interpolations does not change the result, while skipping the interpolation in the analysis only results in large statistical fluctuations.The zero temperature potential for different smearing and interpolation levels for a = 0.0404f m.The label "0× int" means no interpolation used in the analysis.The label "2× int" means that the r-interval used in the interpolation was doubled compared to the default setup.

Appendix C: Analysis of the Wilson line correlators at T > 0
In this appendix we discuss the analysis of the Wilson line correlators at T > 0 Our aim is to gain information on the spectral function corresponding to the Wilson line correlator at T > 0. As discussed in the main text we use the following Ansatz for the spectral function ) where ρ high r (ω) is the dominant part of the spectral function at large ω and is assumed to be temperature independent.Furthermore, ρ peak r (ω, T ) describes the dominant peak in the spectral function and encodes the complex potential at T > 0, while ρ low r is a small contribution to the spectral function below the dominant peak, which is discussed below in more detail.The position and width of the dominant peak in the spectral func-tion should not depend on the interpolating operator details used in the static Q Q correlator, e.g. on the flow time and whether we use Wilson line correlators in Coulomb gauge or Wilson loops.On the other hand ρ low r (ω, T ) and ρ high r (ω) will depend on the specific choices of the interpolating operators used in the correlator, e.g. on the amount of smearing or the gauge tolerance used.In Fig. 9 we show the effective masses for T = 305 MeV and r = 0.606 fm for different flow times.We see non-monotonic behavior and flow time dependence for small τ as we do for the T = 0 case.However, for an intermediate τ -range 0.1 fm < τ < 0.45 fm, where the contribution from ρ peak r (ω, T ) is the dominant one, the effective masses for different flow times agree with each other very well.At τ > 0.5 fm the contribution from ρ low r (ω, T ) becomes important, and we see some dependence on the flow time.As discussed in Ref. [10] ρ low r (ω, T ) depends on the overlap of the chosen Q Q operator with the light states that propagate backward in the Euclidean time together with the forward propagating Q Q.Similar dependence on the level of spatial link smearing of the effective mass was observed in Ref. [10].As discussed in the main text the effective masses corresponding to W sub (τ, r, T ) decrease monotonically with τ , and for sufficiently small τ they are approximately linear in τ .This is demonstrated in Fig. 9, where the effective masses from W sub (τ, r, T ) are shown for T = 305 MeV and r = 0.606 fm at various flow times.Thus the removal of the high energy part of the spectral function also removes the artifacts induced by the gradient flow.For a small contribution from ρ low r (ω, T ) this linear behavior of the effective masses in τ for small τ could be eas- However, the Gaussian form of the spectral function is not physically motivated and the width of the Gaussian cannot be interpreted as ImV (r, T ).
If we assume that the detailed shape of the spectral function away from the peak position is not too important we can define ImV (r, T ) as the width at half maximum height.In this case, a Gaussian form of the spectral function can be used.A physically appealing choice of ρ peak r (ω, T ) is a Lorentzian form.However, this form is only valid for ω values that are not too far from ω = ReV (r, T ).The HTL spectral function of static Q Q [21] is Lorentzian only in the vicinity of the peak and decays exponentially when |ReV − ω| is larger [21].The same holds for the spectral function in the T-matrix approach [18].Therefore, we use a cut Lorentzian for ρ peak r (ω, T ) in our analysis It turns out that the cut Lorentzian also gives an almost linear dependence in τ for the effective masses.In our analysis, we set Cut = 2Γ L .To cross-check our results we also use the Gaussian form.
It was shown in Ref. [10] that the rapid nonlinear decrease in the effective masses is due to ρ low r (ω, T ).This contribution to the spectral function arises from the light states in the medium propagating backward in time which are coupled to the static Q Q propagating forward in time [10].This contribution also depends on the details of the Q Q correlators, e.g.whether one uses Wilson line or Wilson loops and the amount of smearing used [10].We model this part of the spectral function with a single delta function because such a simple form is sufficient to describe the data for the Wilson line correlators with the exception of one data point very close to  the boundary τ = 1/T .We perform fits of subtracted Wilson line correlator with cut Lorentzian form of ρ peak r (ω, T ) a single delta function for ρ low r (ω, T ) for all available data sets omitting the first datum, which is possibly affected by the distortions due to smearing, and the last data point.Some sample fits are shown in Fig. 9 for T = 305 MeV, r = 0.606 fm, β = 7.825, and in Fig. 10 for T = 352 MeV, r = 0.28 fm, β = 8.249.The fits work well as demonstrated in Fig. 10 (bottom), where the relative difference between the fit and the lattice data is shown.Fits using the Gaussian form for ρ peak r (ω, T ) work equally well as demonstrated in Fig. 11.
The amplitude and the position of the small delta function that parametrizes ρ low r are shown relative to the dominant peak in Fig. 12 for β = 8.249 and different temperatures.As one can see from the figure, the position of this delta function is between 1.8 GeV and 3.8 GeV below the position of the dominant peak, and shows only mild dependence on r.The amplitude of this delta function on the other hand increases rapidly with increasing r.Similar results have been obtained for the two other β values.We also note that for small values of r, typically smaller than five times the lattice spacing, it is not necessary to include this small delta function in the fits, i.e. we can set ρ low r to zero and obtain good fits.In Fig. 13 we show the width of the spectral function defined as the width at half of the maximum height as a function of r and different temperatures obtained from the fits using Gaussian and cut Lorentzian form for ρ peak r (ω, T ).We see that using the Gaussian results in a systematically larger width.The Lorentzian parameter Γ L though is dependent on the cut on the Lorentzian.This means that there is a systematic uncertainty in the determination of ImV (r, T ) from the parametrization of the spectral function.As we discuss in the section below it is possible to define the width in a model independent way by considering cumulants of the spectral function.
We also studied the dependence of our results on the real and imaginary part of the potential on the number of smearing level and on the interpolations.As in the zero temperature case we performed the analysis without using interpolation or doubling the interpolation range.We find that the distortions in ReV (r, T ) due to smearing are the same as in the zero temperature potential.This is shown in Fig. 15.This means that while smearing can distort the potential at very short distances, it does not affect the temperature dependence of the real part of the potential.Therefore, we show our results for ReV also at the shortest distances in the main text.At these distance we use slightly different values of c Q shown in Fig. 5 to offset the distortions due to smearing.From Fig. 15 we also show that using interpolation does not introduce a bias in our results for ReV (r, T ).The dependence of ImV (r, T ) on the smearing level and on the interpolations is shown in Fig. 16.As one can see from the figure also here we do not see significant dependence on these.

Appendix D: Cumulants of the spectral function
In this appendix we discuss the cumulants of the spectral functions and their relation to the effective mass of the Wilson line correlators.The cumulants of the spectral functions c n are defined as ) where ⟨...⟩ stands for dωρ r (ω, T ).... Cumulants exist if the spectral function has support in a finite ω range, which is the case for the subtracted spectral function In what follows we will discuss the moments of this spectral function.The cumulants of the spectral function are related to the cumulants of the subtracted Wilson line correlators at τ = 0, m n defined as This can be seen by Taylor expanding the exponential in the spectral representation of the subtracted Wilson line correlator Expanding the exponential in Eq. (D4) and comparing to Eq. (D5) we see that: The first cumulant of the Wilson line correlators is the effective mass.The second cumulant is the slope of the effective mass in τ .We calculated the second cumulant of the subtracted spectral function using the Gaussian form and cut Lorentzian form including and excluding the δ function at small ω.The result of this analysis is shown in Fig. 14.We see that the second cumulant of the spectral function is not sensitive whether we use a Gaussian or cut Lorentzian in   our analysis.Furthermore, the second cumulant does not change much if we include or exclude the contribution from ρ low r (ω, T ) that is the small delta peak.We also see that √ c 2 has a similar dependence on r as the FWHM in Fig. 13 but is somewhat smaller.We have estimated the systematic uncertainties due to the Ansatz for the spectral function by varying the fit range for the cut Lorentzian MeV.The label "0× int" means no interpolation is used in the analysis, while the label "2× int" means that the interpolations range was doubled in the analysis compared to the default setup.Fits done from τ = 2 to 17. Fits with no smearing or no interpolation fits have been cut above r/a ≥ 20 due to large errors.The fit for no smearing is only done up to τ = 16.Error bars are purely statistical.
Ansatz with τmin /a = 2 and τmax /a = {N τ −5, N τ − 4, N τ −3}.The effect of this variation for √ c 2 exceeds the difference upon changing the Ansatz (to Gaussian) or the cut of the Lorentzian.Thus, we used the average of the two most outlying results for cut Lorentzian as central values and the full spread as the systematic error estimate, which we have added in quadrature.For large distances this estimate clearly exceeds the statistical errors.ReV (r, T ) is insensitive to these changes within statistical errors.We can also fit our lattice results on the subtracted Wilson line correlator with the following simple form W approx (τ, r, T ) = exp(m 0 − m 1 τ + m 2 τ 2 /2) (D8) in the range τ /a = 2 − N τ /3, where the effec- The label "0× int" means no interpolation is used in the analysis, while the label "2× int" means that the interpolations range was doubled in the analysis compared to the default setup.Fits done from τ = 2 to 17. Fits with no smearing or no interpolation have been cut above r/a ≥ 20 due to large errors.Fit for no smearing is only done up to τ = 16.Error bars are purely statistical.
tive mass is approximately linear.From this fit, we can then estimate the second cumulant of the spectral function and compare it with the determination of c 2 obtained by integrating the model spectral function based on the cut Lorentzian and the small delta function in almost the entire τ range.This comparison is shown in Fig. 17.We see that the two methods of estimating c 2 are in good agreement.This means that defining ImV (r, T ) in terms of c 2 is model independent and robust.
We also calculated the third cumulant of the spectral function using our fitted spectral function based on the cut Lorentzian form.The result on c 3 , which is the measure of skewness of the spectral function, is shown in Fig. 18.We see that −c 3 is close to zero at small r but then rapidly increases with increasing r.For very small distances, r < 5a ρ low r (ω, T ) was not included in the fit, and therefore, c 3 is exactly zero here.Unfortunately, our lattice results are not precise enough to obtain c 3 using fits with Eq. (D8) extended to higher order polynomials in the exponent.Thus at the present level of accuracy, the short τ behavior of the effective masses can be parametrized solely by m 1 and m 2 .

Figure 1 .
Figure 1.The effective masses at T = 0 and at T ≃ 220 MeV for r ≃ 0.7 fm and a = 0.0280 fm (circles), a = 0.0404 fm (squares) or a = 0.0493 fm (triangles).The green symbols correspond to subtracted data.The lines show the fits discussed in the text.Filled (open) symbols represent ms/m l =20(5).
on N 3 s × N τ lattices for physical strange quark mass, m s and two sets of light (u and d) quark mass, m l = m s /5 and m l = m s /20.The latter corresponds to almost physical pion mass, m π = 161

Figure 2 .
Figure 2. The real part of the potential as a function of r at different temperatures.We show results for a = 0.0280 fm (circles), a = 0.0404 fm (squares) or a = 0.0493 fm (triangles).Open symbols for ms/m l = 5 and filled symbols for ms/m l = 20.

Figure 3 .
Figure 3.The estimate of the imaginary part of the potential from the fit using cut Lorentzian form as a function of r or rT for different temperatures.The three panels focus on different temperature ranges.The circles correspond to a = 0.0280 fm, the squares to a = 0.0404 fm, and the triangles correspond to a = 0.0493 fm.Open symbols for ms/m l = 5 and filled symbols for ms/m l = 20.Error bars include a systematic contribution discussed Appendix D.

Figure 4 .
Figure 4.The effective masses corresponding to the Wilson line correlators at r/a = 15, β = 8.249 obtained for different flow times (top).The effective masses for different flow times after applying the additive shift are discussed in the text (bottom).

Figure 5 .
Figure 5.The additive shifts for different flow times as a function of r/a for β = 8.249.

Figure 6 .
Figure 6.The zero temperature potential for β = 8.249 obtained with different smearing levels, including no smearing.

Figure 7 .
Figure 7. (top) Effective mass for Nτ = 64, Nx = 64, r/a = 20 and τF /a 2 = 0.125 for the raw data, compared to an interpolation fit done around r/a = 20 in a range ±∆r/a = 0.9 with a second order polynomial.(bottom) The correlator as a function of distance r at fixed τ = 15a for the same lattice as the top plot.

Figure 8 .
Figure 8.The zero temperature potential for different smearing and interpolation levels for a = 0.0404f m.The label "0× int" means no interpolation used in the analysis.The label "2× int" means that the r-interval used in the interpolation was doubled compared to the default setup.

Figure 9 .
Figure 9.The effective masses for different flow times at T = 305 MeV, r = 0.606 fm, β = 7.825.The bottom panel shows the effective masses for the subtracted correlator.The lines in the bottom panel show the fits discussed in the text.

Figure 10 .
Figure 10.The effective masses for β = 8.249, T = 352 MeV, r = 0.280 fm and the corresponding fits with the cut Lorentzian plus the delta function for ρ low r shown as a line.The bottom panel shows the relative deviation between the fit and the data with the lines indicating the estimated 1 − σ band of the data.

Figure 11 .
Figure 11.The effective masses for β = 8.249, T = 251 MeV, r = 0.392 fm and the corresponding fits with a Gaussian plus the delta function for ρ low r shown as a line.The bottom panel shows the relative deviation between the fit and the data with the lines indicating the estimated 1 − σ band of the data.

Figure 12 .
Figure 12.The amplitude of the small delta function divided by Ar (left) and the position of the small delta function relative to the position of the dominant peak (right) as a function of r.The results are shown at different temperatures for lattice spacing a = 0.0280 fm (β = 8.249).

Figure 13 .
Figure 13.Width at half the maximum height (FWHM) for the cut Lorentzian fit and Gaussian fit.

Figure 14 .
Figure 14.√ c2 for the Gaussian fit (G), compared to the cut Lorentzian fit (L), or to a Gaussian fit without accounting for the low ω structure (G no delta).

Figure 15 .
Figure 15.Real part of the potential for the cut Lorentzian fit for different smearing and interpolation levels for a = 0.0404f m, Nτ = 20, T = 244 MeV.The label "0× int" means no interpolation is used in the analysis, while the label "2× int" means that the interpolations range was doubled in the analysis compared to the default setup.Fits done from τ = 2 to 17. Fits with no smearing or no interpolation fits have been cut above r/a ≥ 20 due to large errors.The fit for no smearing is only done up to τ = 16.Error bars are purely statistical.

Figure 16 .
Figure 16.Imaginary part of the potential for the cut Lorentzian fit for different smearing and interpolation levels for a = 0.0404f m, Nτ = 20, T = 244 MeV.The label "0× int" means no interpolation is used in the analysis, while the label "2× int" means that the interpolations range was doubled in the analysis compared to the default setup.Fits done from τ = 2 to 17. Fits with no smearing or no interpolation have been cut above r/a ≥ 20 due to large errors.Fit for no smearing is only done up to τ = 16.Error bars are purely statistical.

Figure 17 .
Figure 17.The second cumulant of the subtracted spectral function as a function of r determined from the cut Lorentzian form of the spectral function (circles) and from the second order polynomial fit of the Wilson line correlation function in the τ /a range 2 − Nτ /3 (blue band) for T = 251 MeV and a = 0.0280 fm.

Figure 18 .
Figure 18.(−c3) 1/3 as a function of r in temperature units for lattice spacing a = 0.0280 fm and different temperatures.

Table I .
Parameters for the Nσ = 96, β = 8.249, ams = 0.01011 lattice configurations used.The last column shows the flow time used for each Nτ .In this appendix we discuss the analysis of the Wilson line correlators at zero temperature.For the analysis of the Wilson line correlators, it is Appendix B: Analysis of the Wilson line correlators at T = 0

Table II .
Parameters for Nσ = 64, β = 7.825, ams = 0.0164 lattice configurations.The last column shows the range of flow time in lattice units used in the calculations.The numbers in the square brackets indicate the flow time for which the final results in the paper are presented.
TableIII.Parameters for Nσ = 64, β = 7.596, ams = 0.0202 lattice configurations.The last column shows the range of flow time in lattice units used in the calculations.The numbers in the square brackets indicate the flow time for which the final results in the paper are presented.