A hard thermal loop benchmark for the extraction of the nonperturbative potential
Abstract
The extraction of the finite temperature heavy quark potential from lattice QCD relies on a spectral analysis of the Wilson loop. General arguments tell us that the lowest lying spectral peak encodes, through its position and shape, the real and imaginary part of this complex potential. Here we benchmark this extraction strategy using leading order hardthermal loop (HTL) calculations. I.e. we analytically calculate the Wilson loop and determine the corresponding spectrum. By fitting its lowest lying peak we obtain the real and imaginary part and confirm that the knowledge of the lowest peak alone is sufficient for obtaining the potential. Access to the full spectrum allows an investigation of spectral features that do not contribute to the potential but can pose a challenge to numerical attempts of an analytic continuation from imaginary time data. Differences in these contributions between the Wilson loop and gauge fixed Wilson line correlators are discussed. To better understand the difficulties in a numerical extraction we deploy the Maximum Entropy method with extended search space to HTL correlators in Euclidean time and observe how well the known spectral function and values for the real and imaginary part are reproduced. Possible venues for improvement of the extraction strategy are discussed.
I Introduction
Twenty seven years ago Matsui and Satz Matsui:1986dk () proposed the melting of , i.e. the ground state of the vector channel, as signal for the deconfinement transition in heavyion collisions. The recent success of relativistic heavyion experiments Adare:2006ns (); Tang:2011kr (); Chatrchyan:2012np (); Abelev:2012rv () in observing the relative suppression of charmonium and bottomonium serves as further motivation to develop a first principle description of the phenomena.
In the framework of effective field theories, heavy quarks can be described by nonrelativistic quantum chromodynamics (NRQCD) obtained form QCD by integrating out the hard energy scale, given by the rest mass of the heavy quarks. To describe the bound state of two quarks, one can further integrate out the typical momentum exchange between the bound quarks (see Brambilla:2004jw () and references therein), which leads to potential nonrelativistic QCD (pNRQCD). In this effective field theory the bound state is described by a two point function satisfying a Schrödinger equation.
At zero temperature, the potential between a heavy quark and antiquark is defined from the late time behavior of a Wilson loop and can be directly calculated in Euclideantime lattice simulations or in perturbation theory. At small distances, where perturbation theory converges, both results agree Bazavov:2012ka ().
At high temperature, above the QCD phase transition, one might first expect that the problem becomes simpler as the potential is not confining anymore. Actually, this is not the case since even a proper definition of the potential becomes nontrivial. In fact, the presence of a heat bath is most conveniently incorporated in a Euclidean time framework with finite temporal extend. There the Wilson loop depends on imaginary time and needs to be analytically continued to real time. Only from the large realtime, i.e. behavior, the finite temperature potential can be extracted and happens to be complex Laine:2006ns (); Brambilla:2008cx (). Its imaginary part can be interpreted as Landau damping Beraudo:2007ky () and describes the decaying correlation of the system with its initial state due to scatterings in the plasma.
Along the lines presented in Laine:2006ns (), one can compute the potential in finite temperature perturbation theory. This is a demanding task, as resummations need to be carried out in order to cure infrared divergences. To this day the full result is known only to leading order, whereas a short distance expansion has been calculated to higher order Brambilla:2010vq (); Brambilla:2011mk (). Even if higher orders were available, observing the deconfining transition will remain beyond the reach of perturbation theory.
In Ref. Rothkopf:2011db (), a method was proposed to compute the heavy quark potential nonperturbatively from lattice QCD simulations. Starting from the measurement of the Euclidean Wilson loop on the lattice, its spectral function is reconstructed via the maximum entropy method (MEM). The definition of the potential is based on the peak structure of the Wilson loop spectrum.
Previous numerical evaluations however lead to unexpected results: both the real and imaginary part appear to grow linearly at distances where other quantities, such as the free energies already show significant screening effects. This behavior persisted even at temperatures much larger than the QCD phase transition, where on general grounds, one would expect that the confining potential disappears because of Debye screening Digal:2005ht ().
This problem was solved recently Burnier:2012az () by carefully disentangling the different timescales in the problem. Taking into account the remnants of earlytime nonpotential physics, the lowest lying spectral peak was found to deviate from a naive Lorentzian shape through skewing. Extracted values for real and imaginary part based on this functional form result in a potential that is compatible with Debye screening.
In this paper our aim is twofold: at first we wish to ascertain whether fitting of the lowest lying spectral peak indeed suffices to determine the static heavy quark potential, given the spectral function of the Wilson loop or even the gauge fixed Wilson line correlators. Subsequently it is our goal to better understand the challenges facing a numerical determination of the spectral function by Bayesian analytic continuation. Since in the perturbative approach both Euclidean correlator and spectrum are known, the outcome of the numerical reconstruction can be readily compared.
In section II we review the basics of the method of Ref. Rothkopf:2011db () and its improvement introduced in Burnier:2012az (), which form the basis of the extraction of the potential from lattice simulations. From calculations of the realtime Wilson loop as well as gauge fixed Wilson line correlators in section III we determine and investigate the corresponding spectral functions in section IV. While in section V we apply the peak fitting procedure of Burnier:2012az () to the HTL spectra, section VI scrutinizes how well these spectra can be obtained with the maximum entropy method from the HTL Euclidean correlators. Our conclusion in section VII discusses the limitations of the method and points toward further possible improvements.
Ii Heavy quark potential from Euclidean correlators
The description of the interactions between a pair of heavy quarks and antiquarks at finite temperature in terms of a quantum mechanical potential requires the relevant physics to be well separated from the energy scale of pair creation. In particular
(1) 
needs to be fulfilled^{1}^{1}1See for instance Brambilla:2013dpa () for the discussion of the different limiting cases and their physics, which is satisfied exactly in the static limit (). In that case, the propagation amplitude of an infinitely heavy quark pair can be described by a rectangular temporal Wilson loop where are its temporal and spatial extend. This realtime quantity is defined as the closed contour integral over the matrix valued gauge field along the path of the heavy quarks
(2) 
If the scale hierarchy holds, it is permissible to substitute the field theoretical interactions by an instantaneous potential, so that obeys a Schrödinger type equation
(3) 
At late times, on expect the function to become time independent, so that we may define the heavy quark static potential as
(4) 
Due to the complex weighting factor in Feynman’s path integral, we cannot calculate the realtime Wilson loop using lattice QCD Monte Carlo simulations. Instead we have to rely on an analytic continuation of Euclidean time quantities that are accessible by these numerical methods. In order to connect the heavyquark potential and the Euclidean Wilson Loop one introduces a spectral representation of the realtime quantity,
(5) 
where the time dependence now resides entirely in the integral kernel. Note that the function is not just a Fourier transform but can be shown to be a positive definite spectral function Rothkopf:2009pk ()^{2}^{2}2It is important to distinguish this dependent Wilson loop spectral function from the quarkonium spectral function Burnier:2007qm (); Burnier:2008ia (); Ding:2012iy (); Ding:2012sp () representing the physical quarkonium spectrum.. After analytic continuation one observes that only the integral kernel has changed, whereas the spectral function remains the same
(6) 
Using the Maximum Entropy Method (MEM), a form of Bayesian inference, it is in principle possible, albeit challenging, to invert Eq.(6) and thus to extract the spectral function from . Note that the model independent method of refs. Cuniberti:2001hm (); Burnier:2011jq (); Burnier:2012ts () is not directly applicable as the Wilson loop is not periodic. However a similar method could probably be developed from the general results of refs. Viano1 (); Viano2 (). Once we are in possession of the spectral function we can insert Eq.(5) into Eq.(4), which yields Rothkopf:2009pk ()
(7) 
Direct application of this formula in the case of a numerically reconstructed spectral function is very difficult. It is however possible to determine those structures in the spectral function, which dominate the integral in the infinite time limit.
If we suppose that the time independent potential description holds for all times i.e. in equation (3), an intuitive connection between spectral features and the static potential can be established. In this case equation (3) can be solved and the spectral function turns out to be a simple Breit Wigner peak
(8) 
characterized by its peak position and width .
In general the function however is time dependent at early times and one expects that a wealth of structures, different from the simple Lorentzian example, exists in the spectrum of the Wilson loop at finite temperature. Note that if the potential description is ultimately applicable, the function will become time independent at late times and therefore a corresponding well defined lowest peak must exist. This part of the spectrum encodes all the relevant information on the potential and it alone needs to be reconstructed from the Euclidean correlator.
In Ref. Rothkopf:2011db () it was assumed that the lowest peak is solely described by the late time behavior of the potential and is not affected by the time dependence of the potential at short times. It was shown in Ref. Burnier:2012az () that this is actually not the case. The short time dynamics (nonpotential terms, bound state formation) doesn’t just create additional structures at high frequency but also significantly modifies the shape of the low frequency peak. The most general form of this low peak, derived in Ref. Burnier:2012az (), can be written as
(9)  
Note that this result can also be obtained from pNRQCD where arises from the phase of the singlet normalization factors Brambilla:2004jw ()
In order to calculate the potential from Euclidean correlators we thus need to carry out the following steps:

Calculate the Wilson loop at several separation distances for all possible values along the imaginary time axis .

Use Bayesian inference to extract the most probable spectrum for each value of .
In the following section II we prepare a testing ground for this extraction strategy based on analytic calculations of the realtime and Euclidean Wilson loop in the HTL resummed perturbative approach. Since the analytic continuation can be performed explicitly in HTL, item three of the above list can be tested independently from questions arising from possible inadequacies of the maximum entropy method. The availability of both the spectrum and Euclidean data points on the other hand furthermore allows us to check the degree of success of the MEM itself in the form of a realistic mock data analysis.
Iii Correlators from HTL resummed perturbation theory
iii.1 Wilson loop
In perturbation theory, the Wilson loop is calculated as an expansion in the coupling:
(10) 
starting form . The first nontrivial term () contains only a one gluon exchange and is not enough to describe the correct physics for large Euclidean time . To improve this situation, we resort to the usual ’exponential’ resummation Beraudo:2007ky (), noticing that
(11) 
Thus a better approximation for is
(12) 
as it resums all ’ladder diagrams’ and contains the correct leading order () large behavior.
iii.1.1 Leading order term
We now turn to the calculation of , for which we set the direction along the third spatial axis. In hard thermal loop (HTL) resummed perturbation theory, all diagrams contributing to have one HTL gluon running between the lines of the Wilson loop Laine:2006ns ():
The gluon HTL propagator, written in Euclidean space () and covariant gauge reads:
(14)  
while the HTL selfenergies are given in Appendix A and the projectors take the form:
(15) 
Following Ref. Laine:2006ns (), we rewrite the HTL selfenergies as spectral functions,
(16) 
so that we can perform the sum over analytically:
where we abbreviated the dependence of the second term through the function
(18) 
We can write the spatial vector in spherical coordinates and . In an isotropic plasma, the HTL spectral functions and selfenergies depend on only. Integrating over is trivial and the integral over involves
(19)  
where is the sin integral function. Performing the angular integrals and using gives
(20)  
The first line of equation (III.1.1) is linear in , whereas the next lines are proportional to and therefore symmetric around . We will consider these terms separately in the following:
(21) 
iii.1.2 Part linear in
The part linear in is formally divergent. Using dimensional regularization, the result can be read off from Ref. Laine:2006ns (); the first line of equation (III.1.1) hence gives:
(22)  
In the limit this part yields the real part of the potential:
Note that the result is finite (for ) and the divergence at reflects the behavior of the Coulomb potential.
On the lattice, this term behaves differently^{3}^{3}3The difference with dimensional regularization can be traced back to an infinite constant that is removed in the dimensional regularization procedure.. Roughly speaking, the integral is truncated by the lattice cutoff and thus finite. In this case it is easy to see that it vanishes at , which is expected as a Wilson loop without area is equal to unity. For , it decreases quickly and formally goes to in the limit of an infinite cutoff.
This behavior cannot be canceled by the other terms in equation (III.1.1) as they have a different dependence. It should also not be removed as it encodes the Coulomb part of the potential that we want to obtain. To make a connection to the lattice, we therefore introduce a UV cutoff, mimicking the finite lattice spacing. In this case, performing the integral over the momentum from zero to in equation (22) gives:
where are the and integral function. From the UV regularized version of the correlator we get the following potential,
(24) 
which is plotted in Fig. 7 together with the continuum () potential.
iii.1.3 Symmetric part
We calculate here the symmetric part of the correlator containing the lines 24 of equation (III.1.1). The functions receive a contribution from the cuts of if . For the opposite case they vanish except for a function contribution coming from the pole of . In the following we calculate the contribution from the cuts and poles of the transverse and longitudinal selfenergy separately,
(25) 
As before we introduce a cutoff on the momentum to mimic the effects of the lattice regularization.
Cut contributions
Pole contribution form the longitudinal spectral function
Pole contribution from the transverse spectral function
We proceed in a similar way for the transverse spectral function.
Here, the limit does not exist the integral in equation (III.1.3) is linearly divergent (see Appendix B). Note that such divergences were already observed in Burnier:2009bk (); Berwein:2012mw (), where the Wilson loop of maximal time extend is shown to diverge at next to leading order. The leading order divergence found in Eq. III.1.3 has yet a different nature and consistently vanishes for . In dimensional regularization, it can be shown (see appendix C) to match the cusp divergence Korchemsky:1987wg (); Brandt:1982gz (), which in this case gives Berwein:2012mw ().
Here, we are not interested in trying to renormalize the Wilson loop. It is not needed for our purposes as we aim at a comparison with lattice results, which are also not renormalized. It is however interesting to note that these cusp divergences do not contribute to the potential and only make the Wilson loop heavily suppressed for , hence harder to measure with high accuracy. Removing these divergences in the lattice measurements, without affecting the potential would be of great help to improve the accuracy of the lattice data. One strategy deployed to this end could be the smearing of gluonic links Bazavov:2013zha ().
iii.1.4 Imaginary part of the potential
From the symmetric part, we obtain the imaginary part of the potential,
(29) 
As in the end the infinite time limit will be taken, it is sufficient to consider the low frequency part of the integrals,
Performing the time derivative, using equation (A) and approximating for small as well as the identity
(30) 
we get:
which coincides with the expression obtained in Laine:2006ns (); Beraudo:2007ky (); Brambilla:2008cx ().
iii.1.5 Numerical evaluation
To make close connection to actual lattice data with spatial lattice spacing , we choose to fix the cutoff in our HTL calculations to
(31) 
which naively corresponds to the largest momentum accessible under this finite resolution. Based on a numerical evaluation of the remaining integrals in eq. (21,25III.1.3), we can generate an arbitrary large number of datapoints spanning the imaginary time axis, which carry numerical errors of the order of the machine precision only.
Comparing this ideal HTL Euclidean regularized data to actual measurements from a Monte Carlo simulation in Fig. 1, we find a strong qualitative resemblance. Both graphs exhibit three characteristic features, i.e. a suppression region at small together with an upward trend at , both of which are closely linked to the divergences observed in III.1.3. The datapoints at intermediate are the ones encoding the potential. They exhibit nearly exponential behavior for small separation , where also is small but begin to show noticeable curvature for larger separation distances.
After calculating the realtime values (see Fig. 2) using a similar numerical evaluation of the integrals in (21,25III.1.3), it is possible to obtain the function . As shown in Fig. 2, we can explicitly observe the approach of to a constant value and thus the emergence of a simple exponential behavior of the Wilson loop. Note that in Fig. 2 we show times where the oscillatory behavior is clearly visible while a constant value is actually reached for larger . We refrain from attaching any physical meaning to the length of the swingin period, as it is dominated by the same cusp divergences that lead to the suppression of the Euclidean Wilson loop data points.
iii.2 Gauge fixed Wilson line correlator
Cyclic Wilson line correlators (i.e. color singlet Polyakov loops) fixed to Coulomb gauge have been extensively studied on the lattice, both for the determination of the zero temperature potential as well as in investigations into the freeenergy difference between a medium with and without inserted heavy quarks (see for instance Maezawa:2011aa (); Kaczmarek:2005ui ()). Due to the absence of spatial Wilson lines connecting the temporal links, these quantities offer a significantly better signal to noise ratio than the Wilson loop, especially if the multilevel algorithm Luscher:1996ug () is applied.
Besides the technical question of whether the removal of spatial connectors (or e.g. the application of smearing on spatial links) can lead to an improved lattice observable for the extraction of the potential, it is conceptually of interest to understand whether gauge independent information, such as the potential can be extracted from a gauge dependent quantity such as the Wilson line correlators^{4}^{4}4The crucial difference to potential models is that we do not investigate the single point , but it is the full Euclidean time dependence of the gauge fixed correlator that is used to reveal the values of the potential..
We proceed with the determination of the Euclidean time Wilson line correlator analogously to III.1
(32) 
Calculating in leading order hard thermal loop (HTL) resummed perturbation theory we obtain the expression
(33)  
which contains fewer terms than the Wilson loop of (III.1.1).
iii.2.1 Coulomb gauge
In Coulomb gauge, the HTL Euclidean gluon propagator reads
(34)  
where the self energies are the same as in covariant gauge (see Appendix A). Inserting the propagator into the expression (33) for the Wilson line correlator gives
(35)  
We now rewrite the HTL selfenergies as spectral functions, use the formulas collected in Appendix A to perform the sum over and carry out the angular integrations:
We find that the Coulomb gauge Wilson line correlator features a similar structure as the Wilson loop
(37) 
the symmetric expression however being of much simpler form, depending only on the longitudinal HTL spectral function. At this point we can already anticipate that it is these terms present in both Wilson loop and Wilson line correlator, which contribute to the values of the potential. In particular the cusp singularity connected to the transverse spectral function identified in III.1.3 is absent from the above expression.
iii.2.2 Potential from the Wilson line correlator
As in the case of the Wilson loop, a closed expression for the potential can be obtained using
In the infinite time limit one can make use of
(39) 
which leads us to the same result we encountered for the Wilson loop
(40) 
with the imaginary part given by the integral expression
(41) 
From a practical standpoint this result is encouraging, as it tells us that (to leading order in HTL) the information content regarding the potential encoded in the Coulomb gauge Wilson line correlator is the same as the one found in the Wilson loop. If such a relation persisted into the nonperturbative realm, the absence of cusp divergences and with it the improved signal to noise ratio would make this an ideal observable to reconstruct the potential.
iii.2.3 Numerical evaluation
As for the Wilson loop we wish to compare the Euclidean HTL correlator to actual values measured in quenched lattice QCD Monte Carlo simulations. While the symmetric term is finite, the part linear in still requires a regularization. We deploy the same momentum space cutoff as introduced in III.1.2 and set its value to in the following.
The absence of divergences in the symmetric part of the correlator leads to a significantly different behavior along the imaginary times . As can be seen in the top graph in Fig. 3, where we plot the HTL Wilson line correlator and the first five HTL Wilson loops as comparison. The large suppression at early times as well as the upward trend near are almost absent. Hence most of the datapoints actually carry information on the potential.
Interestingly in the case of the lattice QCD Wilson line correlator, the upward trend is still visible between the last and second to last time step. However contrary to the leading order HTL result, where , the values of these two different correlators on the lattice do not agree at .
iii.2.4 Covariant gauge
The Wilson line correlator can be calculated in the covariant gauge as well. The result depends on the gauge parameter , and contains additional end point divergences Aoyama:1981ev (). These terms however do not contribute in the infinite time limit so that the obtained potential is again the same as in the Wilson loop case (40).
Iv Spectral functions from HTL resummed perturbation theory
iv.1 Spectrum of the Wilson loop
The spectral function can be directly calculated from the realtime correlator via a Fourier transform,
(42)  
We start by analytically investigating the low frequency behavior of this function, as it allows insight into the spectral structures that encode the physics of the heavy quark potential and will be used in its extraction in section V. To benchmark the MEM extraction of the spectrum from Euclidean correlators it is however necessary to compare to the full spectrum, which we will determine from Eq.(42) numerically.
iv.1.1 Analytical estimate for the low energy part of the spectral function
Starting form equation (III.1.1), we introduce the momentum cutoff
(43) 
where the argument of the second exponential function reads
For small frequencies, the main contribution to the spectral function (43) comes form small values of in the above integral. Expanding equation (IV.1.1) around gives:
(45)  
All terms with negative powers of are retained in this expression, as they dominate the integral for late times. Note that the imaginary part of the potential appears as an overall factor in the above expression. Within this approximation, the remaining integrals are carried out analytically and we get:
with and denoting a not near specified normalization constant. From this result, we see that the pole of the spectral function indeed resides at and the width of the peak is closely related to the imaginary part of the potential. The result however is not a Lorentzian, but is precisely of the form (9) derived on general grounds in Burnier:2012az (). Note that the phase related to the skewing of the spectral peak is interestingly also given by the imaginary part of the potential
(47) 
iv.1.2 Full spectral function
We proceed to calculate the full spectral function by integrating numerically equation (42). Applying the discrete Fourier transform to the realtime Wilson loop evaluated on a set of points separated by a , we obtain its values for a wide range of frequencies partly shown in Fig. 4.
As expected from the minute values of at small separation distances, the peak one finds is extremely sharp. However it also becomes clear that the amplitude of the peak is rapidly suppressed as increases. At the same time nonpotential contributions related to the divergent terms in the symmetric part of give rise to a huge background structure spanning a wide range of frequencies.
Note that at a step in the otherwise smooth spectral function is visible. This is a manifestation of the momentum cutoff we introduced to regularize the formally divergent terms. At the same time one can observe that the spectrum continues beyond these frequencies, which is a reminder that the cutoff was not imposed on the HTL gluon spectral functions.
iv.2 Spectrum of the Wilson line correlator in Coulomb gauge
Analogously we can obtain the the spectral function related to the real time Wilson loop correlator
(48)  
At leading order in the HTL resummed expansion, we again have:
(49) 
with
The spectral function can then be calculated analytically close to its peak at small frequency, which yields
Surprisingly at leading order in the HTL approximation we find that the skewing characterized by the quantity is exactly the same as for the Wilson loop. Note that the same result can also be obtained in the covariant gauge.
iv.2.1 Full spectral function
The full spectral functions for the HTL Wilson line correlator are plotted in Fig. 5. One immediately realizes from a comparison with Fig. 4 that even though the peak position, width and skewing are equal to the Wilson loop case, the Coulomb gauge spectral function looks quite different. The first major difference is that the amplitude of the lowest lying peak depends much less on the separation distance , the second is the virtual absence of the background terms populating a large frequency range in the Wilson loop case. Both facts are of course related, since their origin lies in the suppression of the Euclidean Wilson loop correlator induced in the presence of cusp divergences.
V The potential from perturbative spectral functions
Now that we are in possession of the full HTL spectra obtained from both the Wilson loop and the Wilson line correlator in Coulomb gauge, we can test whether the knowledge of the lowest lying spectral features alone suffices to reconstruct the values of the interquark potential in practice. To this end we fit the low region of using the functional form (9) and compare the extracted values with the analytically calculated . We show here the fitting of the Wilson loop spectrum only, since its application to gives exactly the same results (the potential and the skewing are the same). In section VI, where the numerical reconstruction of the spectra from Euclidean time correlator data is concerned, the differences in e.g. the background contributions will however play a major role.