Nonperturbative potential for study of quarkonia in QGP

A thermal potential can be defined to facilitate understanding the behavior of quarkonia in quark-gluon plasma. A nonperturbative evaluation of this potential from lattice QCD is difficult, as it involves real-time corelation function, and has often involved the use of Bayesian analysis, with its associated systematics. In this work we show that using the properties of the static quarkonia thermal correlation functions, one can directly extract a thermal potential for quarkonia from Euclidean Wilson loop data. This leads to a controlled extraction, and allows us to judge the suitability of various model potentials. We also discuss the phenomenology of quarkonia in the gluonic plasma.


I. INTRODUCTION
Quarkonia, mesonic bound states of heavy quark and antiquark, have played a very important role in our understanding of the physics of strong interactions. The experimental signatures of some of these states are distinctive, the most iconic being the dileption peak of the vector quarkonia. In the theoretical side, the heavy quark mass, M Q ≫ Λ QCD , leads to simplifications. The earliest insights about properties of quarkonia states were obtained by treating them as nonrelativistic states bound by a color electric potential. The potential suitable for studies of quarkonia has been calculated in detail using numerical Monte Carlo studies on lattice-regularized QCD; see, e.g., Ref.
[1] for a review. The potential remains an important ingredient in a systematic expansion of quarkonia in 1/M Q [2].
The dilepton peaks of the vector quarkonia, in particular that of the J/ψ, have been extremely important signatures of creation of quark-gluon plasma (QGP) in ultrarelativistic heavy ion collisions (URHIC), following the suggestion three decades ago [3] that the screening of the color charge inside QGP will lead to dissolution of bound states. This was made more quantitative in follow-up studies [4]. The early studies used a perturbative Debye-screened form, which is the free energy of a static QQ pair in perturbative QGP. Here m D is the Debye mass, α = 4 3 g 2 4π , and g is evaluated at a scale determined by the temperature T . Nonperturbatively, the free energy of Q −Q pair in plasma was calculated using lattice QCD [5], which was used as a proxy for an effective finite temperature potential. However, in the early days a proper formalism for potential-based study of quarkonia in QGP was missing. In particular, other thermodynamic quantities can be derived from the free energy, e.g., an "internal energy" for the QQ pair [6]; the use of such quantities have also been explored in the literature [7].
A theoretical formalism for an "effective finite temperature potential", that can be used to study experimentally observed quantities like the dilepton rate, was first provided in ref. [8]. The starting point is a point-split version of the dilepton current, where U is a suitable gauge connection such that V µ is gauge invariant, and the angular brackets denote thermal average. Defining the correlation function the spectral function ρ J (ω, r; T ) is defined from its Fourier transform, ρ J (ω, r; T ) = 1 − e −ω/T dt e iωt C > (t, r).
The dilepton rate is proportional to the spectral function of the point current, ρ(ω; T ) = lim r→0 ρ J (ω, r; T ). Since we are interested in heavy quarks with M Q ≫ T, Λ QCD , Eq. (3) simplifies. Going to the nonrelativistic notation Q = ψ χ where ψ, χ are nonrelativistic fields that annihilate a quark and create an antiquark, respectively, and remembering that since M Q ≫ T , the thermal states do not include Q fields, the leading (M 0 Q ) term in an 1/M Q expansion gives If one has a system where the sole interaction term is a potential V ( r) between the quark and the antiquark, then it is easy to show that, to leading order in 1/M Q , C > (t, r) satisfies [9] i ∂ t − ∇ 2 r M Q C > (t, r) = V ( r) C > (t, r).
In our theory where the QQ are interacting with the thermal medium, we can then define a potential by equating the left hand side of Eq. (6) to V (t, r) C > (t, r) (staying within leading order of 1/M Q ), where the interaction effects are summarized in a time-dependent V (t, r). An effective thermal potential, V T ( r), can then be defined in the large t limit, if the limit exists: V T ( r) = lim t→∞ V (t, r).
The potential V T ( r) can be obtained by going to the static limit, where, modulo renormalization factor, C > (t, r) reduces to a Minkowski-time Wilson loop: Tr P e i t 0 dt1A0(t1, r/2) U (t; r/2, − r/2) P e i 0 t dt2A0(t2,− r/2) U (0; − r/2, r/2) and Eq. (6) reduces to which defines our thermal potential [8,9]. Using V T ( r) to calculate C > (t, r) from Eq. (6) will give the resummation of the leading ladder diagrams. A calculation of V T ( r) in leading order hard thermal loop (HTL) perturbation theory gives [8] (z 2 + 1) 2 1 − sin zx zx (9) and V re T ( r) is given in Eq. (1). In Eq. (9) we have absorved a negative sign in the definition of V im T ( r), so that V im T ( r) takes positive values. V re T ( r) corresponds to the usual physics of Debye screening in medium, such that for sufficiently large screening, the bound states will not form. On the other hand, V im T ( r) clearly leads to a broadening of the spectral function peak. It captures the physics of collision with the thermal particles leading to a decoherence of the QQ wavefunction [10,11]. For the quark and antiquark far apart, r ≫ T , V im T ( r) reaches a finite limit αT giving the damping rate of the individual quarks [9].
It is well-known that the perturbative calculation, Eq. (9), is not suitable at temperatures a few times T c , the deconfinement temperature. The aim of this paper is to make a nonperturbative calculation of an effective thermal potential, using numerical lattice gauge theory techniques. Following the insight of Ref. [8], various authors have tried calculating the thermal potential nonperturbatively. In the next section we will outline our strategy. More details, and some discussion on difference from earlier studies, can be found in Sec. IV.
The potential description, Eq. (6), is of course an approximate description of in-medium quarkonia. First, here the QQ pair is treated as an external probe put in an equilibrium plasma. Then (in the perturbative language) it accounts for a subclass of diagrams. At zero temperature, the justification for this is well-understood. At finite temperature, extra scales come into play, making the picture more complicated. A systematic, effective field theory based study of the interplay of these scales has been made in Ref. [12] in perturbation theory. In the hierarchy of scales one gets the potential Eq. (9), where r B is the radius of the bound state and E B the binding energy. For the temperatures of interest in heavy ion collision experiments, this hierarchy of scales is hardly satisfied. The effective field theory version, however, is perturbative and therefore cannot be directly used for phenomenology.
Instead of going through the nonrelativistic potential route, one could instead try to directly calculate the spectral function by studying the Euclidean J µ J µ correlation function and try to extract the spectral function from it. This has been attempted for charmonia [13] and, using NRQCD, for bottomonia [14,15]. Unfortunately, the extraction of spectral function from the Euclidean correlator is a notoriously difficult problem, and the systematics are large (see [16] for a discussion, and [17] for early comparison of potential model results with results of [13]). Therefore a nonperturbatively determined potential continues to be important for quarkonia phenomenology; see, e.g., [18]. In recent years, there have also been attempts to come out of the picture of external probe in equilibriated plasma, by treating the quarkonia in plasma as an open quantum system [11,[19][20][21]. The potential remains an important structure in such frameworks [19,21].
The plan of the rest of the paper is as follows. After explaining the calculational methodology in the next section, in Sec. III we will give the calculational details. Sec. IV will give our results for the potential. Some phenomenological discussions and implications of the potential obtained will be discussed in Sec. V, and the last section will have a summary and discussion.

II. NONPERTURBATIVE STUDY OF FINITE TEMPERATURE POTENTIAL
The potential V T ( r) is directly related to the Minkowski space Wilson loop, Eq. (8). But in numerical Monte Carlo studies we work in Euclidean space. At zero temperature, it is straightforward to calculate the QQ potential from the Euclidean Wilson loop: At finite temperature, the simple spectral decomposition outlined in Eq. (10) does not work. The first attempt to extract the QQ potential from W T (τ, r) was carried out in Ref. [22]. The spectral decomposition of the Minkowski-time loop leads to [22] Bayesian techniques were used to extract ρ(ω, r; T ) from W T (τ, r), and then calculate the potential using Eq. (11). The reconstruction of ρ(ω, r; T ) from W T (τ, r) is a notoriously unstable problem. To make matters worse, the quality of the Wilson loop data deteriorates quickly at large τ (this problem can be somewhat alleviated with recent numerical techniques [23]). While very impressive technological improvements have occurred in the Bayesian analysis techniques, the results obtained for potential still have stability issues or have large errorbars, especially for V im T ( r). The first calculations [22] employed a Bayesian analysis method similar to Maximum entropy and fitted the spectral function peak with a Lorentzian form. The results obtained, however, are substantially different from a later analysis [24] which is of similar philosophy but employs a slightly different Bayesian analysis, and fits to a skew-Lorentzian form [25]. The state-of-the-art for calculations in the gluonic plasma follow a similar methodology and can be seen in Ref. [26]. Studies have also been carried out for full QGP (i.e. with thermal quarks), both with a Lorentzian form of the spectral function [27] and using Bayesian reconstruction methods [28]. While the improvement in the analysis method has been impressive, the results still suffer from stability issues; in particular, it is not easy to disentangle the effects of V im T ( r) and V re T ( r) in W T (τ, r). In this paper we take a different approach. Let us motivate it by writing The physics of V re T ( r) is very similar to that of the zero temperature potential, Eq. (10). We therefore expect the real part of the potential to come from the part of w(τ, r) which has a linear behavior around β/2:w(τ, r) ∼ −(τ − β/2)V re T ( r) + .... We isolate thew part by splitting W T (τ, r) as follows: We find that B(τ, r) = log W a T (τ, r) has exactly the behavior we were expecting: B(τ, r) ∼ β 2 − τ V re T ( r) over a large range of τ around β/2. We illustrate this in Figure 1, where B(τ, r)/(β/2 − τ ) is plotted. We also checked that for configurations below T c , where we can extract the potential from the full wilson loop, W a T (τ, r) gives the same result but reaches the plateau sooner.
In order to understand the behavior of W p T (τ, r), we write a spectral decomposition for A(τ, r) = log W p T (τ, r): To go to the potential, we follow the usual route of going to real time τ → it: The potential is obtained in the large time limit of Eq. (15), when the oscillating factors exp(±iωt) ensure that only the ω → 0 contribution to the integral survives. In this limit exp(βω) → 1 and it is obvious from Eq. (15) that A(it) leads to an imaginary potential. One can then extract the real and imaginary parts of the potential from W a T (τ, r) and W p T (τ, r) respectively [8], The argument above is motivated by perturbative studies of the potential, where the split Eq. (16) has been noted [8]. Even with Eq. (16), it is not obvious that the extraction of the potential from the Euclidean correlation function is simple; Eq. (16) involves large Minkowski time, while the nonperturbative data that can be obtained from the lattice is in Euclidean time τ ∈ [0, β). Successful extraction of potential from Eq. (16) is contingent upon the contribution from the "potential modes" dominating the behavior of the correlation functions A(τ, r), B(τ, r). Fortunately, this is what was found in the behavior of the nonperturbative data. As we already discussed above and showed in Figure 1, over a large range of τ , W a T (τ, r) ∼ exp(−cτ ), leading to a straightforward extraction of V re T ( r) from the slope of the exponent. We actually obtained very similar plateau in all our lattices. See Sec. IV A for more discussion.
One, of course, does not expect such a simple behavior from A(τ, r): Eq. (14) rules out a simple linear behavior near β/2. This is expected: if W p T (τ, r) had a linear exponential falloff, it would have contributed to a real potential! The large time behavior of W p T (τ → it, r) can be inferred from a closer examination of Eq. (15), using the fact that in the limit of large t, exp(−i ω t) − exp(i ω t − ω β) −→ −2π i ω δ(ω). Then in order to get a potential Interestingly, this leading singularity structure gives a very good qualitative description of the τ dependence of ∂ τ A(τ, r). This is illustrated in the right panel of Figure 1. The argument in this section is based on the assumption that a thermal potential can be defined using Eq. (8). We then make plausibility arguments on the structure of W T (τ, r), and show that the nonperturbative lattice data supports this structure. The arguments leading to Eq. (16) can be made more concrete using Feynman diagrammatic language [9]: in Appendix B we outline this argument. There we also show the results of the leading order HTL perturbation calculation of W T (τ, r) [8], which fully supports the structures of A(τ, r) and B(τ, r) discussed above, and which motivated this nonperturbative study. The 1/ω 2 behavior in Eq. (17) comes from the term ρ(ω) ω 2 and a distribution function, (1 + n B (ω)) ω→0 − −− → T ω , which follows from the structure of the time-ordered propagator (see Appendix B and Eq. (B4)). It is connected to the scattering origin of the imaginary part of the potential, discussed below Eq. (9).
Our strategy for extraction of the potential is therefore straightforward: we extract V re T ( r) from a linear fit to B(τ, r) and to get V im T ( r), we expand σ(ω; T ) in Eq. (14) in the basis (1 + n B (ω)) {1/ω, ω, ...}, and extract V im T ( r) from the coefficient of the most singular term. As Figure 1 suggests, the leading terms dominate the data around β/2, allowing us to extract the potential relatively simply. We discuss further details in Sec. IV.

III. TECHNICAL DETAILS OF OUR STUDY
In this work, we have calculated the QQ potential in a gluonic plasma, for moderately high temperatures ≤ 2T c . We have generated lattices with a space-time anisotropic discretization with ξ = a s /a τ ≈ 3. A convenient algorithm for doing this is given in [29]. We follow this reference to estimate the lattice parameters we require. The anisotropy is estimated nonperturbatively from comparison of spatial and temporal Wilson loops [29], while a τ is estimated from the string tension calculated from temporal Wilson loops. We use three sets of lattices, with a τ ranging between 1/19T c and 1/45T c . For each set, we change the temperature by changing N τ , while keeping the spatial volume fixed.
For each set, we first make short Monte Carlo runs at closely spaced N τ to find the N τ for deconfinement transition. The final lattice sets used for the studies above T c are shown in Table I. For much of this paper, we will measure all scales in units of T c . However, for Sec. V we will need to quote physical units. We will do so by taking the string tension σ = (0.44 GeV) 2 . This translates to a transition temperature ∼ 280 MeV. The spatial extent of the lattices are 1.44 fm or above. Some more details regarding the runs are given in Appendix A.
In order to determine the potential, we calculate thermal expectation values of timelike Wilson loops, i.e., the Euclidean time version of W M in Eq. (8). It is well-known that for the spatial connections U straight thin-link gauge connections are not suitable: they lead to very noisy signals in numerical Monte Carlo studies. To alleviate the problem due to extended spatial connections, we do APE smearing [30]. This constitutes of a replacement of the elementary gauge links U i , iteratively. The spatial connections U are then constructed from these smeared links. For this work, we have taken α = 2.5. Note that Eq. (18) does not involve the time direction, and the time direction links are not smeared. So time slices and the definition of transfer matrix is not affected by the smearing. We use the multilevel algorithm [23] in the temporal direction: this allows us to get a good signal even for Wilson loops with large time extent. For calculation of the potential at T = 0, smearing is routinely used, and the potential should be independent of the smearing. In the finite temperature case, the extracted "potential" may depend on the details of the connection U; but the actual physical quantity one is interested in, the quarkonia peak in dilepton channel, is independent of it, as it is connected to the point current. We do, however, do a detailed study of the dependence of the potential on the smearing level in the next section.
In the literature, the correlator of Coulomb gauge fixed Wilson lines have often been used to extract the potential. The Coulomb gauge fixing can be formally understood as a dressing of the quark fields [31]: where ψ Ω (x) = Ω(x)ψ(x) and Ω(x) is a dressing function such that ψ Ω (x) is gauge invariant [31]. The Coulomb gauge potential has obvious advantages in that the extended spatial links are not there. At T = 0, it is also easy to argue (and has been well-tested) that the Coulomb gauge potential agrees with the potential extracted from the Wilson loop. For T > T c such detailed comparison does not exist in the literature. Here we have made such a comparative study. The coulomb gauge is fixed to an accuracy of 10 −7 . We have also checked that the results  do not change if the accuracy is made 10 −6 or 10 −9 instead. The potential from this Wilson line correlator has also been presented in Sec. IV. In particular for the imaginary part of the potential, we observe differences between this potential and that obtained from the smeared Wilson loop. Since the Wilson loop operator does not involve dressing of the quark field, the connection to the point-point correlator at r → 0 is transparent. We use the potential obtained from the Wilson loop for further studies in Sec. V.

IV. POTENTIAL CALCULATED FROM WILSON LOOPS
In this section we present the details of our extraction of the potential, using Eq. (16). In Sec. IV A we discuss the real part of the potential. The results for the free energy of a QQ pair is given in Sec. IV B, and the extraction of V im T ( r) is discussed in Sec. IV C. Besides quoting the results for the potential, we also compare the potential at different levels of smearing, and the results for Coulomb gauge. Finally, in Sec. IV D we will discuss the spectral representation Eq. (11), and touch on issues of direct extraction of spectral function from Euclidean data.
A. Real part of the potential As outlined in Sec. II and Figure 1, for smeared Wilson loops the extraction of the real part of the potential from W a T (τ, r) is straightforward. Defining a local potential through −∂ τ log W a T (τ, r) shows a plateau near β/2. In the left panel of Figure 2 we show the "local measurements" of V re T ( r) from Wilson loops at different levels of smearing. The errorbars shown are from a Jackknife analysis, after blocking the data to reduce autocorrelation. As the figure shows, while for a small number of APE smearing steps, the local mass takes time to reach a plateau, on increasing the number of steps a plateau is reached quickly, and we can easily extract the potential using a single exponent fit. While we have shown the local mass for one particular case, the effects are very similar for all our sets. For each smearing level the value obtained from the fit is shown by the horizontal band of the same color. The goodness of the fit, as demonstrated by χ 2 , is very good. The figure also shows that varying the number of smearing steps over a large range does not seem to have any statistically significant effect on the value reached at the plateau.
We also show in the figure the local values of the potential obtained from the Coulomb gauge fixed Wilson lines. As the figure shows, the Coulomb gauge data seems to be noisier than the data from Wilson loops. We checked that this is not an artifact of the accuracy at which the Coulomb gauge is fixed. Also the Coulomb gauge results are found to be close to the results from the smeared Wilson loops, but the difference between them is statistically significant.
In the right panel of Figure 2 we summarize the fitted value of V re T ( r) for this set. At this scale, the dependence of the potential on the smearing level is hardly visible. Similarly, the potential from smeared Wilson loops and that from Coulomb gauge fixed Wilson lines are very close, though they differ at 1σ level.
As we have discussed in Sec. III, we believe that for study of quarkonia property in medium, the potential from the smeared Wilson loop is appropriate. It is satisfactory that V re T ( r) becomes practically independent of the level of smearing very soon. Anyway, when quoting a result for V re T ( r), we include, as a systematic error, some variation with the level of smearing: for example, for the set shown in Figure 2 we include the spread in results between smearing  levels of 100 and 250 as a systematic error. In what follows, our error bars for V re T ( r) include this variation for all sets.
Results from lattices at a finite lattice spacing have discretization errors. We can have an idea of the size of the discretization error by comparing the results at different lattice spacings. As Table I shows, we have lattices with three different lattice spacings at 1.2 T c , and at 1.5 T c we have results with two different lattice spacings. In Figure  3 we show the potentials calculated from lattices at different lattice spacings. Within our error bars the results agree very well, indicating that the cutoff effects are very small at these lattice spacings. We will, therefore, take the results on our finest lattice spacings as a valid estimator of the continuum results. Figure 4 summarises our results for V re T ( r) at different temperatures. We see that the potentials at the two temperatures below T c agree completely, indicating that the temperature effect is small even at 0.75 T c . The potentials have the familiar Cornell form, with a dip at small r and a linearly rising part for r 0.5 fm. This behavior changes abruptly on crossing T c : while the short distance part, 0.2 fm, remains similar to the form below T c , beyond rT c ∼ 0.5 ∼ 0.35 fm the effect of string breaking clearly shows up, and the potential becomes flatter with increasing temperature.

B. Free energy
The study of the free energy cost of introducing a QQ pair in the plasma is almost as old as the study of deconfinement transition in QCD. The free energy of QQ pair was calculated from the correlator of Polyakov loops, L( r) L † ( 0)  The inset highlights the long distance part. The Coulomb gauge result is seen to be close to that obtained from smeared Wilson loop, but with statistically significant difference. [32]. Later, the free energy cost of a singlet QQ pair was connected to the cyclic Wilson loop (for sufficiently smeared loops) [33]: or from Coulomb gauge fixed Circular Wilson lines [34] (see also [35]). In leading order perturbation theory, the singlet free energy agrees with V re T ( r). The singlet free energy has been studied in great detail, for both gluonic plasma and the theory with quarks [5], and we do not intend to add to the existing results. Here we will, however, examine the issue of whether the perturbative agreement between the free energy and V re T ( r) is also valid nonperturbatively. In Figure 5 we show the singlet free energy calculated from the smeared circular Wilson loop at different levels of smearing, and that from the Coulomb gauge fixed operator. The smearing dependence is similar to what was seen for V re T ( r): the results are quite insensitive to the smearing level used. The Coulomb gauge operator is close to the Wilson loop results, but not exactly identical.
In Figure 6 we compare the free energy and V re T ( r), extracted from the smeared Wilson loops, at three different temperatures. As discussed before in Sec. IV A, the results are expected to be valid continuum results. At all temperatures, we find that F ( r; T ) and V re T ( r) are very close to each other. However, at long distances V re T ( r) shows slightly less screened behavior than F ( r; T ). C. Imaginary part of the potential As we have discussed in Sec. II, the behavior of the symmetrized correlation function W p T (τ, r) is dominated by the most singular behavior in Eq. (17), which is the term that corresponds to V im T ( r). Encouraged by this, we expand σ(ω; T ) in Eq. (14) in a series where the form of Eq. (21) is motivated by the structure of A(τ, r) (see Eq. (B4) and the discussion at the end of Sec. II). σ(ω; T ) in Eq. (21) has the property that σ(−ω; T ) = e −βω σ(ω; T ) and so the integrand in Eq. (14) is an even function of ω; the even powers of ω are absent in Eq. (21) as they won't contribute to the integral. The imaginary potential V im T ( r) is obtained from the coefficient of 1/ω term: V im T ( r) = π β c 0 . Putting Eq. (21) in Eq. (15), we get the linear series for the "local mass": where the generalized ζ functions ζ(s, x) = ∞ n=1 1 (x + n) s . Note that this form Eq. (22) is similar to, and could also be motivated by, perturbation theory [8].
The data near β/2 gives a very good fit to just two terms in Eq. (22), and with three terms, almost the entire range of τ could be fit in all our data sets. In Figure 7 we show the results for V im T ( r) obtained with different levels of smearing. The error bar here includes the variation due to change in number of terms of Eq. (22) in the fit. The dependence on the level of smearing is stronger here, but a plateau can be reached after some levels of smearing. When quoting a result for the imaginary part of the potential in what follows, our error bar encompasses the spread among the different smearing levels in this plateau.
In Figure 8 we show the imaginary potential at two different temperatures, obtained on lattices with different cutoffs. While our coarsest lattice, set I, seems to show some lattice spacing dependence, the results from the two finer sets agree very well. We therefore take V im T ( r) obtained from our finest lattice as a good approximation to the continuum result.
In Figure 9 we show our final results for the imaginary potential at three different temperatures. In Sec. V we will use this data as the nonperturbatively evaluated V im T ( r), and explore its physics. We have shown here the results above T c only; we have, however, run the same analysis strategy on the configurations below T c , and checked that the results are consistent with zero, as expected.  Combining the results of Sec. IV A and Sec. IV C, we can write the correlation function near the center of the lattice as where the higher order terms, do not contribute to the potential. For explaining the Wilson loop data over a substantial range near the center, just c 1 is enough, while adding c 2 allows us to explain W T (τ, r) over the entire range except a couple of points at the edge. Further insight into the potential can be obtained if we investigate the structure of the low ω part of ρ(ω, r; T ) in Eq. (11). In order to do this, we take the Fourier transform of the structure of W T (τ, r), Eq. (23), continued to real time:W T (t = −iτ, r). This shows a peak structure at low ω, as has been anticipated in various lattice extractions of the potential, e.g., [22,[24][25][26]28]. Interestingly, however, the peak structure is very different from what has been often anticipated. In the literature often a Lorentzian or a Gaussian structure has been assumed for the peak. Instead, we find a structure that is exponentially falling in the low ω side of the peak, ∼ exp(ω/T ), while in the high ω side it falls only like a power law. Illustration of the peak structure is shown for a few representative values of r in Figure  10. Given this peak structure, we could rephrase our discussion of the potential extraction by simply starting from a structure like those shown in Figure 10, and extracting the potential from them. We checked numerically that the laplace transform of the peak gives a statistically satisfactory description of W T (τ, r) near β/2. While the direct Bayesian inversions have to grapple with the issue of convergence of the integral in the negative ω side, here we could easily do the integral by putting a lower cutoff: because of the sharp fall, the effect of the cutoff on the value of the integral is negligible. The addition of the correction terms do not have any significant effect on the position or the half-width of the peak, but modifies the fall-off with ω away from the position of the peak. Bayesian statistics based studies of the potential proceed without making strong assumptions about the structure of the peak. In fact, some of the Baysian analyses use only very mild information about the peak. We would like to add a note of caution here. If we do not make the assumption Eq. (22), which is well-motivated by the physics involved in the imaginary potential and also by perturbation theory, it is possible to describe the Wilson loop data by other structures, leading to different V im T ( r). In particular, a very good description of the data is provided by the form The spectral peak obtained from this form is considerably different from that shown above; see Figure 11. A Bayesian analysis, in our opinion, ought to include the broad features of the low ω peak discussed in the previous paragraph.

V. DISCUSSION OF POTENTIALS AND QUARKONIA
Let us try to analyze in some detail the potentials obtained in Sec. IV. We start with V re T ( r). Figure 4 shows our estimation of V re T ( r) at different temperatures. As is expected for a gluonic plasma, the thermal effects are negligible at temperatures of 0.75 T c : the potential agrees completely between 0.75 T c and 0.63 T c . So the potential at our lowest temperature measured for each set can safely be taken to approximate the zero-temperature potential. The potential shows the familiar features of the 1/r singularity at short distances and the linear rise at large distances, and gives a good fit to the Cornell form.
As we cross T c , the finite temperature potential is close to that at T =0 at short distances. But clear temperature effects are seen as r increases: in particular, the linear behavior of the T = 0 potential gets screened. In perturbation theory one expects, in leading order, a Debye-screened form of V re T ( r) which is same as the free energy [8], where m D = gT in leading order and α(T ) is the running coupling at the appropriate temperature scale. In Figure  12 this form, Eq. (25), is shown at different temperatures, along with the nonperturbatively obtained potential. For drawing the perturbative curve, following [8], we have used one-loop formula for the coupling [36], α −1 (T ) = 33 8π log(6.742 T /Λ MS ), and T c /Λ MS = 1.10-1.20 [37]. The band in the perturbative form in Figure 12 corresponds to this range in T c /Λ MS . Since we are interested in the r dependence of V re T ( r), the additive renormalization constant C is fixed by matching to the lattice potential at rT c = 0.5 at T = 2T c .
As Figure 12 shows, the perturbative form does not explain the potential obtained in Sec. IV A. In particular, the long distance part of the potential is not as flat as the screened Debye form predicts: as if a shadow of the string tension rise survives.
Since the long distance part of the QQ potential in QCD vacuum has a linear string tension term, a natural next step would be to try a screened form of the string tension term. The string tension term being entirely nonperturbative, there is, however, no single unique/preferred possibility for the screened form of this term. We will consider here two models for screening that have been discussed in the literature. A linear string tension is obtained in the 1+1 dimensional Schwinger model. Since string is essentially an one-dimensional object, one can assume that the physics of screening of the string term will also be similar to that in the Schwinger model. Such a consideration leads to the   (31).The errors shown include the variation with fit range, and should be treated as a systematic band rather than a statistical 1 − σ band.
This form of the screened potential can also be obtained by generalizing the timelike gluon propagator to [39] The second term gives a linear string term in the limit m D → 0. We treat Eq. (27) here as a purely phenomenological construct to model the screening in string tension term. In perturbation theory, one expects α to be a function of r and T . In the Cornell potential, however, one usually treats α as a constant. We follow [38] and keep α, σ fixed to their T = 0 value, the temperature dependence entering in Eq. (26) only through m D . The long distance part of the potential V re T ( r), Figure 4, is fitted to Eq. (26) to obtain m D , C ′ . V re 1D ( r, T ) does a good job of explaining the measured potential as shown in Figure 12. We have tried a few fit ranges covering the large distance side of our measured potential. The band in Figure 12 shows the variation of the fit parameters on shifting the fit range. The narrowness of the band is evidence for the stability of the fit to the form of Eq. (26). The fitted value of m D obtained from the fits is shown in Table II; the range corresponds to this change in fit range.
A different line of argument to a screened potential is to start with a generalized Gauss' law which gives a linear potential [40]. The medium effect then can be incorporated by introducing a medium permittivity [41]. Using an isotropic permittivity motivated by HTL perturbation theory leads to the potential [42] V re 3D ( r, T ) = − α r e −mDr − Γ(1/4) 2π where µ 2 = m D σ α , x = µr and K 1/4 is the modified Bessel function of the second kind [43]. At large r, the second term behaves like exp −x 2 /2 √ x . The results of the fit to this form are also shown in Figure 12 and the value of m D shown in Table II. The fit to Eq. (26) is found to be slightly more stable than that to Eq. (28), and so we use it for analysis of quarkonia behavior. However, Eq. (28) also approximately captures the r dependence of V re T ( r); with our data we can not statistically rule out either of the one-dimensional and three-dimensional screening forms.
The imaginary part, V im T ( r), turns out to be more difficult to model using the conventional screening forms available in the literature. The perturbative form of the imaginary part, Eq. (9), is shown in Figure 13 together with our data, for three different temperatures. The parameters used are identical to that for the real part, as detailed below Eq. (25). The data shows very different behavior from that of Eq. (9): at short distance, the perturbative result overshoots the data, but it soon saturates, while our nonperturbative data does not show a sign of saturation in the distance scale studied by us. The perturbative result V im T ( r) pert behaves ∼ r 2 log r at small r, and saturates to ∼ α T as r → ∞. The nonperturbative data shows a r 2 behavior to a much larger distance: in particular, almost the whole range of r explored by us, rT c 1, can be fitted to a quadratic behavior at 1.5 T c and 2 T c .
The HTL permittivity that leads to Eq. (9), is complex, so as to produce a complex potential. Use of this permittivity in the generalized Gauss' law leads to [42] Eq. (28) and an imaginary part V im 3D ( r, T ) is also shown in Figure 13, with the legend '3D'. Here the value of m D obtained from Eq. (28) is used, and the band corresponds to the range in m D (Table II). This form has a similar behaviour ∼ r 2 at small r to the data. While it is steeper at large r than the perturbative form, it is less steep than our data.
If one uses a complex permittivity analogous to the HTL term in conjunction with the modified propagator of Eq. (27), one can get the "complex potential" for 1D screening, i.e., the imaginary part of Eq. (26). The imaginary part reads [39] V im 1D ( r, T ) is shown in Figure 13 with legend '1D'; the value of m D is that obtained from Eq. (26) in Table II. This form seems to have a higher slope than our data at small r and a smaller slope at large r, though at 1.5 T c it is close to our data in the range of r studied by us.
As Figure 13 reveals, none of the simple forms discussed does a good job of modelling our data for the imaginary potential over the range of r studied by us. At small r, the numerically calculated potential has a smaller slope than either the screened string forms or the forms Eq. (30) and Eq. (29). At large r, on the other hand, it is steeper. Both of these latter forms, in turn, show a much larger imaginary part than the perturbative form at large r, with Eq. (29) comparable to our data at larger values of r.
We are interested in the ground state quarkonium peaks in the spectral function. While it is most sensitive to the short distance part of the potential, it is also affected by the long distance part, especially as the binding energy becomes less and the state becomes broader. As we mentioned before, in the range rT c 1 studied here, our data for V im T ( r) grows ∼ r 2 . Of course, on physical principles we expect it to saturate at large r. Motivated by Eq. (30), we tried to model the imaginary part of the potential by fitting the data to an arbitrary combination of V im pert ( r, T ) and V im σ ( r, T ). We also tried to fit it to a purely quadratic form. Finally, we calculate the spectral function for finite mass quark through integrating Eq. (6), with where V re 1D ( r, T ) is given in Eq. (26) and the parameters b, a 1 , a 2 are given in Table II. We emphasize that our forms for V im fit ( r, T ) represent purely phenomenological fits of the data; one can take them to correspond to two limiting asymptotic behaviors given the data. We will treat the results for the spectral function obtained with the two forms of V im T ( r) in Eq. (31) as a systematic band, and look for features of the band. In the left panel of Figure 14 we have shown the spectral function obtained this way at 1.5 T c , with the quark mass varying from 1.5 GeV to 6 GeV. At T = 0, using the unscreened Cornell potential we get a series of sharp peaks. We denote the mass of the 1S state as M P , and normalize the x axis with respect to it in Figure 14. At 1.5 T c , even for a quark mass of 6 GeV we only find one peak. Of course, 1.5 T c here corresponds to a temperature of about 420 MeV. Expectedly, the peak is the sharpest for the heaviest quark, gradually broadening till, for quark masses close to the charm, only a very broad peak structure can be seen. The spectral function for M Q = 1.5 GeV is also qualitatively different from the others, and is very different from the spectral function obtained directly from J/ψ correlators in [13], but in qualitative agreement with a later study [44]. Similar results have been seen in [17]. In the right panel of Figure 14 we have shown the results for the Υ peak. Quark mass M Q was tuned to get the 1S meson mass ∼ 9.45 In the x axis, the zero is at ω = MP , the 1S state mass obtained using the Cornell potential.
GeV. A sharp peak is seen at 1.2 T c , which gradually broadens as the temperature increases. But a peak structure survives all the way to 2 T c . Note that 2 T c here corresponds to about 560 MeV, setting the scale using the string tension. Also at low temperatures the peak is quite narrow, in comparison to what was found from nonrelativistic bottomonia correlators in [14].

VI. SUMMARY
One popular way of studying the medium modification of quarkonia in quark-gluon plasma is through defining an effective "in-medium" potential. Theoretically, a suitable potential can be defined [8,9] by examining the time dependence of the Minkowski-time Wilson loop, Eq. (8). This potential is complex, with the real part of the potential describing the Debye-screened binding of the QQ pair in medium and the imaginary part related to damping of the wavefunction due to interaction with the thermal medium. A nonperturbative extraction of this potential [22] involves extracting the low-frequency structure of the spectral function from the Euclidean-time Wilson loop, Eq. (11). This is in general a very difficult problem. The existing studies in the literature have either progressed through using Bayesian analysis, with their associated, and sometimes hard-to-estimate, systematic errors, or by making ad-hoc assumption about the low-frequency structure.
In this work we have introduced a new method of nonperturbative evaluation of the potential. We find that a reorganization of the Euclidean Wilson loop data, motivated by the underlying structure of the finite temperature correlation function (see Appendix B), leads to an enormous simplification in the extraction of the peak structure from the Wilson loop. The main ingredients of our method are outlined in Sec. II, and the details are given in Sec. IV. We have calculated the potential in a gluonic plasma for temperatures ≤ 2T c from smeared Wilson loops, calculated using anisotropic lattice discretization of the gluonic theory. Our results for the potential are summarized in Figure 4 and Figure 9. The real part of the potential, which shows the Cornell form below T c with no noticable temperature dependence upto 0.75 T c , shows Debye screening on crossing T c , with the screening mass increasing with temperature. The form of the potential is different from the perturbative form at least upto 2 T c , as illustrated in Figure 12. The imaginary part of the potential is zero below T c . Above T c it rises rapidly, with a spatial dependence ∼ r 2 till distances r 1/T c . Its behaviour is sharply different from the perturbative result, as illustrated in Figure  13.
In the course of our study, we have also investigated issues like dependence of the finite-temperature potential on the definition of the operator, which, we feel, have not been properly discussed in the literature. We have examined how the potential depends on the smearing, and compared the potential obtained from smeared Wilson loops with those from Coulomb gauge fixed Wilson line correlators. We have also examined the relation between the real part of the potential and the free energy of a static QQ pair in the plasma; see Figure 6. In Sec. IV D we have discussed the structure of the low energy peak of the spectral function. It is quite different from the Lorentzian structure that has often been assumed in direct extractions of potential from Euclidean wilson loop using Eq. (11). We have also illustrated, with example, the difficulty of extracting the low energy peak from the Euclidean Wilson loop without putting in additional physics input.
Our data for the extracted potential can be found in Sec. IV, in particular in Figure 4 and Figure 9. Moreover, for various purposes it is convenient to have a parametrization of the potential. In Sec. V we have explored various standard forms of a screened potential. As Figure 12 shows, for the real part, the form of 1D screening of the string potential seems to give a reasonable description of the data, with parameters given in Table II. For the imaginary part it is more difficult to find quantitative agreement with a standard screened form. The potential rises ∼ r 2 till intermediate distances rT c ∼ 1. While the potential is expected to saturate as r → ∞, it is difficult to make any statement about that behavior from our data at rT c 1. A purely phenomenological generalization of Eq. (30), using arbitrary linear combination of V im T ( r) pert and V im T ( r) σ , seems to give a good description of the data in the range of r explored by us, with coefficients given in Table II. Since we expect the long distance behavior of V im T ( r) to be somewhere between this and the r 2 behavior, for a study of quarkonia in the plasma we use both of the forms for V im T ( r) (see Eq. (31)). The difference in the spectral structure obtained with these two forms is considered as a qualitative systematic band.
The spectral function peaks for S-state quarkonia with different quark masses are shown in Figure 14. In the left panel, the variation of the spectral function with M Q is shown. Below T c the spectral function has a number of narrow peaks corresponding to the nS states; but above T c only the 1S peak survived even for M Q = 6 GeV. For M Q = 1.5 GeV, close to the mass of charm, there is no significant peak structure at this temperature. Of course, the nonrelativistic formalism may not be valid for charmonia at these temperatures. For M Q = 3 GeV, a clear peak structure is seen at 1.5 T c , with very little shift in the peak position. In the right panel of Figure 14 the spectral function for 1S bottomonia is shown. While the peak structure weakens with temperature, a clear peak survives till 2 T c , with very little shift in peak position, and reasonably narrow peak, at least till 1.5 T c .
While potential by itself does not provide a complete description of medium interaction of quarkonia, it is an important part of a complete description, and can provide useful inputs for more sophisticated studies like direct extraction of spectral functions; they also can provide essential nonperturbative ingredients of an open quantum system analysis of in-medium quarkonia [19,21]. Our results are for the quenched theory, and one needs to be careful when applying them for quarkonia phenomenology. They, however, provide benchmarks for comparing with direct extractions of quarkonia spectral functions from euclidean correlators. More importantly, the method we have outlined for the extraction of the potential, Sec. II, is quite simple and stable, and we expect one should be able to use it to extract reliable potential also from dynamical lattices.