The imaginary part of the heavy-quark potential from real-time Yang-Mills dynamics

We extract the imaginary part of the heavy-quark potential using classical-statistical simulations of real-time Yang-Mills dynamics in classical thermal equilibrium. The r-dependence of the imaginary part of the potential is extracted by measuring the temporal decay of Wilson loops of spatial length r. We compare our results to continuum expressions obtained using hard thermal loop theory and to semi-analytic lattice perturbation theory calculations using the hard classical loop formalism. We find that, when plotted as a function of mDr, where mD is the hard classical loop Debye mass, the imaginary part of the heavy-quark potential shows little sensitivity to the lattice spacing at small mDr ≲ 1 and agrees well with the semi-analytic hard classical loop result. For large quark-antiquark separations, we quantify the magnitude of the non-perturbative long-range corrections to the imaginary part of the heavy-quark potential. We present our results for a wide range of temperatures, lattice spacings, and lattice volumes. This work sets the stage for extracting the imaginary part of the heavy-quark potential in an expanding non-equilibrium Yang Mills plasma.


Introduction
At extreme energy densities nuclear matter undergoes a phase transition from a state characterized by confined hadrons to a state in which quarks and gluons become the relevant degrees of freedom. Ultrarelativistic heavy-ion collision experiments at the Relativistic Heavy Ion Collider in New York and the Large Hadron Collider at CERN have now collected a wealth of data concerning the properties of the high energy density phase of nuclear matter, which is called the quark-gluon plasma (QGP) [1]. One of the key observables used in measuring the properties of the QGP is the suppression of heavy quark-antiquark bound states such as bottomonium. The suppression observed by experimentalists gives a measurement of the in-medium breakup rate of bottomonium states and allows one to constrain key plasma observables such as the initial central temperature of the QGP generated in heavy-ion collisions, its shear viscosity to entropy density ratio, and the differential rate at which the QGP expands in non-central collisions (see e.g. [2][3][4][5]). Fundamentally, the computation of the survival probability of a given bottomonium state can be cast into the framework of open quantum systems (OQS) in which there is a probe (bottomonium states) and medium (light quarks and gluons). Within the OQS framework, in order to describe the in-medium evolution of bottomonium states one must trace over the medium degrees of freedom and obtain evolution equations for the reduced density matrix of the system [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. In the limit that the medium relaxation time scale and the intrinsic time scale of the probe are much smaller than the probe relaxation time scale, the resulting dynamical equation for the reduced density matrix can be cast into a so-called Lindblad form [20,21]. A key outcome of such calculations is that the effective heavy-quark potential possesses an imaginary part which can be related to the total inmedium decay width of the states. This imaginary part has been determined using direct quantum field theoretic or effective field theory calculations [10,12,22].
There have been computations of the imaginary part of the heavy-quark potential based on high-temperature quantum chromodynamics (QCD) calculations in the hard thermal loop (HTL) limit [22][23][24][25][26][27][28][29], using effective field theory (pNRQCD) [30][31][32][33], finite-temperature lattice QCD [34][35][36][37][38][39][40][41][42][43], and real-time classical-statistical solutions of Yang-Mills theory in classical thermal equilibrium [44,45]. In this work we build upon the studies presented in Ref. [44] and present findings that are complementary to Ref. [45]. In Ref. [44] the authors presented first results for the imaginary part of the heavy-quark potential using classicalstatistical Yang-Mills simulations on spatially 3D lattices of size 12 3 and 16 3 . In this paper we extend these results to larger lattices up to 252 3 and consider SU(2) and SU(3) gauge theories. Additionally, the results of Ref. [44] were presented only for a few spatial points r in a table. Due to the use of rather large lattice sizes, we can now compute the imaginary part of the heavy quark potential at larger values of r/a and reconstruct the functional form of the imaginary part of the heavy-quark potential for a much wider range of distances. This allows us to make more precise comparisons between our lattice-extracted imaginary part and (a) results obtained in the continuum limit using hard thermal loops [22] and (b) results obtained using the lattice-regularized hard-classical-loop (HCL) theory [44].
Herein we will present results for the imaginary part of the heavy-quark potential obtained using classical Yang Mills (CYM) simulations of a thermalized gluonic plasma. The use of CYM simulations is motivated by the fact that in situations where (a) gluonic occupation numbers are large, such as in thermal equilibrium for sufficiently low momenta or in the initial stages of heavy-ion collisions, and (b) the gauge coupling is weak g 2 1, vacuum contributions to observables are suppressed by powers of the gauge coupling. As a consequence, neglecting vacuum contributions in such cases often leads to a good approximation to the full quantum dynamics.
However, one should note that in the real-time classical lattice theory rotational invariance is broken and the short and long distance physics do not decouple [46,47]. As a consequence, the results that we obtain in our classical simulations may differ from those in thermal equilibrium.
In order to extract information about thermal systems using CYM, one must prepare the thermalized field configurations. Historically, this is done by preparing 3+1D configurations using Monte-Carlo techniques [48][49][50][51][52][53]. Herein we follow a simpler strategy and initialize the system close to thermal equilibrium using momentum-space initialization and let the fields thermalize in real time before extracting observables. The advantage of this procedure is that we can run simulations on very large lattices with small lattice spacings a at moderate computational cost. In practice, however, this means that we have to extract the temperature g 2 T numerically in the thermalized system from correlation functions. We find that this only introduces a small uncertainty for g 2 T of the order of a few percent, which justifies this approach.
One additional complication that arises when dealing with CYM treatments is that they do not possess a finite ultraviolet limit due to the Rayleigh-Jeans divergence and the theory is non-renormalizable [52,[54][55][56][57][58]. For this reason, it is important to identify a suitable manner in which to scale results in order to extract relevant information. We will use the Debye mass m D computed in the hard-classical-loops framework to demonstrate that, when plotted as a function of m D r, the imaginary part of the heavy-quark potential is only mildly sensitive to the lattice spacing, or more generally to the simulation parameter β ∝ 1/(g 2 T a), at small distances m D r 1. Additionally we find that, when presented in this manner, results obtained using CYM simulations agree well with the semi-analytic hard classical loop result at small quark-antiquark separations. Our study also suggests that the latter approaches a finite form in the large-β limit for all separations. We present our results for a wide range of temperatures, lattice spacings, and lattice volumes, which map to values of β in the range 16 → 300.
This work sets the stage for extracting the out-of-equilibrium imaginary part of the heavy-quark potential in expanding Yang Mills plasmas [59,60]. Classical Yang-Mills simulations are applicable in that case, allowing one to extrapolate to small lattice spacings, because the dynamics of the highly occupied plasmas are governed by hard scales that can be set below the lattice momentum cutoff. Such non-equilibrium extractions are necessary since it is currently unknown how to analytically extract the heavy-quark potential in anisotropic plasmas due the presence of a non-Abelian plasma instability called the chromo-Weibel instability [25,[61][62][63][64]. In particular, recent studies have shown that particularly in highly anisotropic systems nonperturbative effects beyond hard loop calculations are crucial [65,66].
The structure of our paper is as follows. In Sec. 2 we provide the theoretical background necessary for the computation of the imaginary part of the heavy-quark potential. In Sec. 3 we present the results from lattice simulations and the comparison to perturbative calculations. We present concluding remarks and outlook for future works in Sec. 4.

Theory and numerical setup
We consider pure SU(N c ) gauge theory with the Yang-Mills classical action with Einstein sum convention for repeated Lorentz indices µ, ν = 0, . . . , 3. The field strength tensor is given by with gauge coupling g, gauge field A µ (x) and commutator [, ]. Unless stated otherwise, we will use N c = 3.

Lattice discretization and equations of motion
We use a standard real-time lattice discretization approach where fields are discretized on cubic lattices with N 3 sites and lattice spacing a (see, e.g., Refs. [60,67] and references therein for more details). In this real-time approach, spatial gauge fields are replaced by gauge links U j (t, x) ≈ exp (igaA j (t, x)) at discrete coordinates x k = n k a for n k = 0, . . . , N − 1, while temporal gauge with A 0 = 0, and thus U 0 = 1 is used. The classical equations of motion are written in a gauge-covariant manner as and are solved alternately in a leapfrog scheme [68]. The plaquette is defined as defines the anti-Hermitian traceless part of the matrix V . The time step is taken to be small, typically a t /a = 0.01 in this work, to reduce temporal lattice artifacts. The classical equations (2.2) guarantee that the degree of Gauss law violation is conserved at every time step Moreover, this lattice discretization guarantees gauge invariance of observables like Wilson loops that are supposed to be gauge invariant.

Computation of observables
Classical observables are computed by averaging over different configurations evolved independently and their uncertainity corresponds accordingly to the uncertainty of the mean.

Imaginary part of the classical potential
We are interested in extracting the imaginary part of the classical potential V cl (t, r) with r ≡ |x|. Following Refs. [22,44], it can be calculated using as the asymptotic temporal slope of log[C cl (t, r)]. The classical thermal Wilson loop C cl (t, r) is defined as with temporal Wilson lines W [(t 0 , x); (t, x)] = 1 and spatial Wilson lines W [(t, 0); (t, x)] = U j (t, 0)U j (t, a j )U j (t, 2 a j ) · · · U j (t, x) for x =â j r andâ j = a j /a being a spatial unit vector.
Since the classical thermal state is homogeneous, the Wilson loop is additionally averaged over all lattice points by averaging over the reference coordinates 0.
In order to extract the imaginary part of V cl (t, r) we compute the time-dependence of C cl (t, r). Due to the imaginary part of the in-medium heavy quark potential, this quantity will decay exponentially at late times with the rate of exponential decay set by the imaginary part of V cl (t, r). As a result, V cl (t, r) can be computed using the logarithmic slope of this decay as mentioned above. We note that this is formally applicable only in the large-t limit, but we find that an exponential decay is established rather quickly, particularly when using large lattice sizes. We then define the imaginary part of the static classical potential as the late-time limit (2.7) In practice, this agrees with the logarithmic slope extracted from C cl (t, r).

Construction of quasi-thermalized and fully-thermalized classical states
In order to associate our measurements of V cl (t, r) with a specific temperature T , it is necessary to have a method for self-consistently determining the temperature of the CYM fields. One can construct thermalized CYM configurations using 3+1D Monte-Carlo techniques [48][49][50][51][52][53]. However, as mentioned in the introduction, we use a simpler technique, which amounts to initializing the fields in a quasi-thermal configuration in momentum-space, as will be explained around (2.10), and then allowing them to self-thermalize dynamically. On large lattices this method is quite efficient and one finds that the fields self-thermalize quickly. One complication is that the resulting equilibrated temperature T is different than the quasi-thermal initial temperature T 0 used to initialize the fields in momentum-space. As a result, one must have a method to extract the time-dependent temperature in order to (a) determine when field thermalization is functionally complete and (b) extract the final thermalized temperature of the CYM fields. Once the system is fully thermalized (within an acceptable uncertainty), one can then proceed with the measurement of V cl (t, r). For this purpose, we extract chromo-electric field correlation functions in Fourier space. In practice, the chromo-electric fields are Fourier transformed as and projected onto normalized transverse vectors v λ j (p) with two polarizations λ = 1, 2 transverse to the longitudinal polarization vector v 3 j (p) = p j /p, providing E λ (t, p). Here we used x k = n k a with k = 1, 2, 3, the lattice momentum definition p k a = −i(1 − exp(−2πim k /N )), p 2 = k |p k | 2 , and n k , m k = 1, . . . , N . The correlation functions are then computed as and similarly for A fields. This expression implies averaging over the direction of p due to the system's approximate isotropy for not too large momenta and we will consider EE T /L as functions of p. 1 Note that these correlation functions are, in general, gauge dependent. We get rid of the residual gauge freedom by gauge-fixing to the Coulomb-type gauge ∂ j A j (t, x) = 0 to almost machine precision right before computing the correlators. As mentioned above, while all measurements will be performed in classical thermal equilibrium, for practical purposes we start our simulations in a state very close to thermal equilibrium and let the system dynamically thermalize. We use two-point correlation functions EE in Fourier space to check to what extent the system has thermalized. Our procedure is as follows.
We set the fields in Fourier space at initial time t = 0 such that they satisfy These initial conditions can be implemented efficiently by setting (see, e.g., [60,67,69]) with the adjoint color index a = 1, . . . , N 2 c − 1, and with Gaussian distributed complex random numbers satisfying (α λ a (p) while all other correlations vanish. Subsequently, we restore the Gauss law D j E j (x) = 0 with covariant derivative D µ to machine precision using the algorithm in Ref. [51]. This algorithm changes the chromo-electric field correlators mostly at low momenta and, for instance, makes EE L finite. Then we let the system evolve according to the Yang-Mills classical equations of motion (2.2). Classical thermal equilibrium is reached when the EE T /L (t, p) correlation functions do not visibly change over an ex- An example of the results obtained using this thermalization algorithm is shown in the left panel of Fig. 1 for EE T and simulation parameters g 2 T 0 = 0.375, N = 64, a = 1, averaged over N config = 224 configurations. The dashed orange horizontal line shows the initial correlator prior to Gauss law restoration, that equals g 2 T 0 for all momenta. After the restoration, the transverse correlator shows some momentum dependence and is given by the blue dashed line at t/a = 0. In the subsequent evolution, the system thermalizes rather quickly. Already for times t/a 30, the correlation function is almost stationary, for later times t/a 60, the deviation between curves at different times lies within the uncertainty of the data, as can be seen in the inset. At that point, the system can be considered approximately thermal.
As shown in Ref. [70] and recently used in [71], the correlation functions EE T /L (t, p) encode the temperature and the Debye mass of the classical state (2.12) These relations are not exact but become more reliable at larger momenta p m D . We can use them to extract the temperature g 2 T and an estimate for the Debye (screening) masŝ m D . A systematic error (for g 2 T ) is estimated as the deviation of the correlators from these simple functional forms, as will be explained shortly. The extraction of these parameters from the transverse and longitudinal chromo-electric field correlations is demonstrated in the right panel of Fig. 1 for the same parameters as in the left panel at the time t/a = 140. One can see that the thermalized correlations are lower than the initial value for a = 1 (left) and a = 0.1 (right). All parameters are shown in lattice units. g 2 T 0 for the considered parameters and that they indeed agree well with the fit functions , as shown by the black dashed and purple continuous lines. Remarkably, these independent fits lead to almost the same value of the temperature and we set Since the transverse correlator EE T still grows slowly with p, we extract the value g 2 T UV at the largest momentum to estimate the error as ∆T = T UV − T . This systematic error is shown by the shaded band in the figure. We find that the uncertainty in the extraction of the temperature is generally quite low and becomes smaller with decreasing lattice spacing. For the coarse lattice spacing a = 1 used in Fig. 1, the relative error is only ∆T /T < 3%. The extracted temperature and mass values, g 2T EtEt , g 2T E l E l andm D , are shown in Fig. 2 as functions of g 2 T 0 for the larger lattice spacing a = 1 (left) and a fine lattice spacing a = 0.1 (right). One observes that both g 2T EtEt and g 2T E l E l , visible as blue and red points, respectively, agree well in both cases. While for the small lattice spacing, g 2 T 0 and the extracted g 2 T are almost identical, one observes deviations that grow with temperature for the coarser lattice. Since the lattice simulations depend on the combination β = 2N c /(g 2 T a), as will be explained in Sec. 2.6, we can summarize that deviations increase with decreasing β, while for growing β values, both ∆T /T and |T 0 − T |/T become smaller, enabling a better control of the value of the temperature.

The HTL, HCL, and lattice Debye masses
We will also need to compute the Debye mass in order to compare results obtained using different analytic and numerical approaches. The definition of the Debye mass depends on the method being used. In this paper we consider three methods for defining the Debye mass: (1) continuum hard-thermal-loops calculations, (2) lattice-discretized hard classical loops calculations, and (3) direct lattice measurement using chromo-electric field correlators. The definitions in these three cases are:
2. In classical thermal equilibrium one has f BE (p) → T /p, and the discretization on an anisotropic cubic lattice with a finite lattice spacing leads to the following hard classical loop (HCL) expression where Σ ≈ 3.1759 is a factor that results from the anisotropy of the lattice [52]. Strictly speaking, this expression is the leading-order perturbative result in the infinite volume limit. Since we employ large lattice sizes, the latter effect can be neglected. Moreover, finite volume effects are expected to be less important than sub-leading perturbative effects that typically contribute to (2.14) and (2.15) at low momenta p m D . We will use (2.15) in our figures when we plot the CYM lattice results as a function of m D r.
3. Our final possibility is to measure the Debye massm D directly from our thermalized CYM configurations using chromo-electric field correlators, which is explained in the previous subsection. Its temperature, and thus β dependence is shown in Fig. 2 for different lattice spacings. One generally finds that it grows with temperature. More precisely, we find an approximately linear connection to m HCL for 92 ≤ β ≤ 300, while at lower β deviations from a linear connection become more sizable.

Parameters and scaling properties with the coupling
The coupling g can be scaled out of the dynamical equations of motion (2.2) by rescaling all field amplitudes as gA → A and gE → E. As visible in (2.12), the coupling then only enters in the combination g 2 T . This dimensionful scale can then be used to make all quantities dimensionless by rescaling them with appropriate powers. In the rescaled version, one then has EE T → 1. Note that this rescaling implies that lattice simulations only depend on the lattice size N 3 and on the lattice spacings g 2 T a ≡ 2N c /β and g 2 T a t . Since we use a t a, temporal artefacts that may result from a dependence on g 2 T a t are suppressed, and the simulations mainly depend on If not stated otherwise, we will write our values for a, g 2 T 0 , the extracted temperature g 2 T and all dimensionful variables in lattice units. We will also provide the corresponding β values.

Results
In this section we present the non-perturbative results from classical-statistical simulations of real-time Yang-Mills theory on cubic lattices with N 3 sites and the spacing length of a.
As detailed in Sec. 2, the fields are initialized with initial temperature T 0 , they thermalize by solving classical equations of motion, and the actual temperature T is extracted with an error estimate from correlation functions of chromo-electric fields. We present the values of the imaginary part of the classical potential Im[V cl (r)] for SU (3) extracted from the evolution of temporal Wilson loops [22,44]. For comparison, we start this section by recalling the corresponding analytic results from dimensionally regularized HTL and the expressions from lattice-regularized HCL calculations of Im[V cl (r)] to the second-order. We compute the HCL potential numerically for a wide range of β values and fit the numerical HCL results to the functional form of a suitable analytically available approximation with only two parameters. The fit allows us to extract the β dependence of the parameters of the HCL potential, including an estimate for the large-β limit with β → ∞. We then discuss our numerical results from classical Yang-Mills simulations for a wide range of β values, and compare them with previously published data of Im[V cl (r)] obtained from lattice calculations, with simulations in SU (2) theory and with HCL results. Our main simulation results are also summarized in tables in App. A.

Hard Thermal Loop result
In Ref. [22] the authors derived an expression for the imaginary part of the heavy-quark potential to leading-order in the strong coupling constant using the continuum hard thermal loop framework and dimensional regularization. Their final result could be expressed compactly as and m HTL D given by the continuum hard thermal loop Deybe mass (2.14). Note that in the large-r limit one has lim r→∞ φ(r) = 1. As a result the asymptotic value of the imaginary part in this case is −C F g 2 T /4π.

Hard Classical Loop results
The result for the imaginary part of the classical potential in the infinite volume and infinite time limit (N → ∞, t → ∞) from second-order perturbation theory regularized on a cubic  lattice (HCL: Hard Classical Loop) of size (aN ) 3 is calculated in Ref. [44] and is given by where the accented variables are defined as We perform the numerical integration necessary in (3.3) using the VEGAS algorithm [72], as implemented in the CUBA library [73]. The results for g 2 T = 0. 44 Fig. 3 as functions of m D r with m D given by the HCL Debye mass (2.15). As can be seen from this panel, the asymptotic large-r value of Im[V HCL ] depends on β, however, the small-r behavior is similar for all values of β shown. We find that the numerical HCL curves are fit very well by the functional form To make the agreement better visible, we show in the right panel of Fig. 3 the same results divided by their corresponding asymptotic values A HCL ∞ . All curves fall on top of each other, demonstrating that Im[V (2) cl (r)]/(g 2 T A HCL ∞ ) is well described by φ B HCL m HCL D r . This also shows that, while fitted values of A HCL ∞ depend on β, this functional form and, in   Since the large-β limit corresponds to the limit a → 0, its existence is very nontrivial for an observable in classical thermal lattice field theory due to the Rayleigh-Jeans divergence. Our analysis suggests that such a limit exists for Im[V regularization of Sec. 3.1. As mentioned before, these are The A ∞ values are quite close, with the estimated HCL value |A HCL ∞ (β→∞)| being 35% larger than the corresponding HTL value. Note here that our extrapolation is expected to have a large uncertainty for A HCL ∞ for β → ∞ since its value is still quite rapidly changing in the considered β interval, as seen in Fig. 4. While the fit works surprisingly well over a large β interval, which gave us confidence to write a limiting value, our analysis is not a proof that |A HCL ∞ (β→∞)| is indeed finite but rather an observation from this fit function. This is different for the parameter B, that agrees well with its estimated large-β limit for β 200. Since B multiplies the divergent Debye mass, the product B m HCL D could be interpreted as a replacement for the lattice regularized mass.

Nonperturbative CYM simulation results for the classical Wilson loop
We now turn to the results of our CYM lattice simulations. To extract the imaginary part of the classical potential Im[V cl (r)] from non-perturbative lattice simulations, we follow the procedure introduced in Ref. [44] and revisited in Sec. 2.3. We consider rectangular Wilson loops of different spatial size r = na and temporal length t = n t a t . 3 We use a temporal lattice spacing of a t = a/100 for all results reported in this paper. The imaginary part of the classical potential Im[V cl (r)] at separation r for a heavy quarkonium system is calculated from the large-time slope of the logarithm of the Wilson loop C cl (r, t), averaged over lattice sites, orientations and configurations. The time evolution of C cl (r, t) is shown for two example cases in Fig. 5.
In the left panel of Fig. 5 we present typical results obtained on a small lattice which is averaged over N config = 28 classical-statistical configurations. The bands indicate the uncertainty in the mean obtained by dividing the standard deviation of the mean value by N config N iso N start.points . As indicated from this last expression, we exploit the isotropy and homogeneity of the lattice in order to obtain increased statistics, where we average over two different Wilson loop orientations and use N start.points = N 3 possible spatial starting points from which we can measure the Wilson loop. As can be seen from this panel, one cannot perform the fits to extract the logarithmic slope at early times. In practice, one must wait a sufficient amount of time beyond the thermalization time (t 0 ). Typically, we use a fit window that encompasses (t−t 0 )/a > 10. We find that as long as one uses starting times for the fit window which satisfy this constraint, there is only a very small variation in the extracted logarithmic slopes.
In the right panel of Fig. 5 we present typical results obtained on a large lattice in which we only take into account one configuration for the initial fields. As this panel demonstrates, when using large lattices, the uncertainty is decreased due to the much larger number of starting points that can be used to measure the Wilson loop. As a result, on large lattices one does not need to sample as many initial configurations, which helps to reduce the run time required to extract statistically accurate results. We have explicitly checked that the uncertainty bands obtained with one initial configuration are consistent with results obtained using larger number of initial configurations but a smaller lattice size. In Fig. 6 we compare our results to the previously available ones from Ref. [44] for β = 16. We set the parameter g 2 T 0 = 0.44 for the initial distribution so that the resulting temperature extracted from EE correlators is g 2 T ≈ 0.373, which is close to the value    Fig. 6 include uncertainties from averaging the Wilson loops over the starting points, orientations, classical configurations, uncertainties coming from the fitting residuals, and uncertainties from the temperature extraction. As can be seen from this Figure, our results are in agreement within uncertainties with the reported results from Ref. [44].

Comparing results for SU(3) vs. SU(2)
To check for the N c dependence of the SU (N c ) gauge theories, we perform here similar lattice calculations for SU (3) and SU (2) systems and extract Im[V cl (r)]/g 2 T for different values of T 0 . In Fig. 7 we use the same lattice parameters a = 1 and N = 12 as we did in Fig. 6 and compare the simulations for the values g 2 T 0 = {0.2, 0.3, 0.4}. In the left panel of Fig. 7, the values of the ratio Im[V cl (r)]/(g 2 T C F N c ) are shown for the two gauge groups as functions of r/a. One finds that they agree well at low r/a while the deviations grow with the distance. However, the agreement improves for smaller temperatures (larger β values).
In the right panel, the ratio of the imaginary parts of the potentials of the SU (3) and , is plotted as a function of r/a. The horizontal dashed line corresponds to (N c C F [SU (3)])/(N c C F [SU (2)]) = 8/3, which is the ratio of their N c C F values. One observes that the data points are close to this value for all r and T 0 , which shows that the imaginary part of the potential in our lattice simulations admits Casimir scaling This is similar to the perturbative HCL expression (3.3) where due to β ∝ N c , one also finds this relation (3.9). The deviations visible in the right panel of Fig. 7 decrease with increasing β (smaller temperature) and with decreasing distance, as we have observed in the left panel. An important source for these deviations is that we compare here simulations with the same g 2 T 0 parameter and not with the same temperature g 2 T . As visible in the left panel of Fig. 2 for a = 1 and SU (3) theory, the extracted temperatures deviate stronger at larger g 2 T 0 . The deviations in SU (2) theory are qualitatively similar but can lead to quantitatively different values of g 2 T , which implies that we compare the SU (N c ) theories effectively at slightly different temperatures. This effect is reduced when going to smaller temperatures or smaller lattice spacings (larger β values). Moreover, as we will show below, the short-distance behavior is less sensitive to the temperature, which is the reason why the deviations decrease at lower r/a. As can be seen from this Figure (Fig. 8), the CYM lattice results seem to approach a finite large-β form when plotted as a function of m D r, with the CYM results with β = 136 even overlapping with the corresponding HCL curve at m D r 1. At large values of r our CYM simulation results approach the corresponding HCL curve for the respective β value. For the largest β = 136 in the figure, we know from the left panel of Fig. 4 that the amplitude of the potential at large r is still quite far away from the extrapolated value. However, Fig. 8 indicates that with increasing β, the CYM potential is approaching the perturbative HCL large-β limit from below.
We emphasize that, while previous studies showed the potential for different β values as functions of r/a, we find that plotting it as a function of m D r makes the comparison more intuitive. Most importantly, this allows to study the classical potential even at large β values while incorporating the dominant UV divergence into the mass. 4 We have performed such a large-β extrapolation for the HCL calculations already in Sec is an estimate of the Debye mass including higher-order and nonperturbative contributions. We find that, in practice, the HCL Debye mass (3.3) is proportional to the correlatorextracted Debye massm D over a wide range of β via Eq. (2.16) with proportionality constant 1.42, and that the resulting rescaling only introduces a constant rescaling of the horizontal axis. Since there are larger uncertainties inherent in the extraction ofm D from the chromo-electric correlators, we have chosen to rescale the classical lattice and HCL results using the same HCL mass definition.

Small-distance behavior
Let us now compare the fitting function (3.5) with parameters A ∞ and B to the imaginary part of the heavy-quark potential. We saw previously that it agrees with this form in both the HTL and the HCL formalisms with the parameters We provide evidence that, at small distances, the classical potential extracted from our lattice simulation data also follows the functional form (3.5) for different values of β. This is shown in the left panel of Fig. 9, where Im[V cl (r)]/g 2 T for fixed lattice spacing a = 0.1 is plotted for different temperatures corresponding to 100 β ≤ 300, as a function of m D r. Fits to each data set using (3.5) in the considered interval m D r ≤ 6 are included as continuous lines. One finds very good agreement between the data and the respective fits. Moreover, although m D and the extracted potential vary with the temperature T and, equivalently, with β, all curves are seen to fall on top of each other at low distances m D r 1. We have seen such behavior already in Fig. 8 for even smaller values of β. 5 This indicates approximately temperature-independent low-distance behavior for a wide β interval. Averaging over the resulting parameter values from the fits to each data set in Fig. 9, we obtain A latt ∞ = −0.24 ± 0.02 = (−0.18 ± 0.015) C F , B latt = 1.2 ± 0.05 . (3.12) Here the errors have been estimated by adding the standard error and the fitting error. In the right panel of Fig. 9, we combine all of the curves of the left panel and compare them to the function (3.5) with the mean parameter values in (3.12), shown as a black line. As expected, one observes very good agreement with our data, especially at low distances m D r. We can also perform a short-distance expansion of the fitting function for the imaginary part of the potential (3.5), neglecting terms O((m D r) 4 ), The resulting curve is shown as a red line and is observed to agree well with our data points for m D r 0.5. Thus, for a wide range of β values, the short-distance behavior of our CYM lattice data agrees well with the perturbative functional form (3.5) and its leading short-distance expansion, which is parametrically given by |Im[V cl (r)]| ∼ C F g 2 T (m D r) 2 log(m D r).

Conclusions and outlook
In this paper we used classical-statistical lattice simulations of pure Yang-Mills fields to extract the imaginary part of the heavy-quark potential as a function of the quark-antiquark separation. In order to carry out our simulations on large lattices, we used a simplified scheme to generate thermalized gauge field configurations which relied on initialization of chromo-electric fields in momentum-space followed by a period of self-thermalization. In order to determine whether the fields were thermalized, we extracted the temperature of the system as a function of time directly from the dynamically generated chromo-electric field correlators.
We went beyond previous classical-statistical lattice calculations of Im[V cl ] by considering rather large lattice sizes and a wide range of values of β (or, equivalently, of the lattice spacing a and temperature). We considered both SU (2) and SU (3) gauge groups and found that the potentials obtained obeyed a simple Casimir scaling. Using our CYM simulations on such large lattices, we argued that Im[V cl ] should be plotted as a function of m D r, with the screening mass m D , when comparing our lattice results for different β values or to perturbative calculations. In particular, we showed that Im[V cl ] approaches the perturbative lattice-regularized HCL results with increasing β. As a side result, our calculations suggest that a finite large-β limit exists for Im[V cl ] when calculated in the HCL framework and plotted as a function of m D r. In the region of small distances m D r 1, we demonstrated that Im[V cl ] from our CYM lattice simulations is insensitive to the finite values of β over a wide parameter range, and is close to the HCL results.
We also found that both our lattice simulation and HCL results were very well approximated by a functional form which can be obtained from a leading-order hard-thermal loop calculation. One only needs to take into account a different prefactor A ∞ and a different scaling of the argument of the functional form, encoded in a parameter B. Using fits of this form and then expanding the result at small m D r, we were able to extract small-distance approximations for the imaginary part of the heavy-quark potential.
Looking to the future, we plan to perform a similar calculation in an anisotropic gluonic plasma that is expanding along the longitudinal direction. This case is of phenomenological interest since the expansion of the quark-gluon plasma in relativistic heavy-ion collisions generates large momentum-space anisotropy. This renders the perturbative analytic calculation of the imaginary part of the heavy-quark potential ill-defined due to the presence of unstable modes in the soft gauge field propagator. In contrast, classical-statistical lattice simulations are applicable in that case, and this work sets the stage to corresponding future studies out of equilibrium.  Table 2. Values of Im[V cl (r)]/g 2 T for different T 0 extracted from lattice simulations with a = 0.1, N = 252 (one configuration).