Proper static potential in classical lattice gauge theory at finite T

We compute the proper real-time interaction potential between a static quark and antiquark in classical lattice gauge theory at finite temperature. Our central result is the determination of the screened real-part of this potential, and we reconfirm the presence of an imaginary part. The real part is intimately related to the back-reaction of the static sources onto the gauge fields, incorporated via Gauss's law. Differences in the treatment of static sources in quantum and classical lattice gauge theory are discussed.


Introduction
The bound states of heavy quarks and antiquarks, so called heavy quarkonium, constitute a unique laboratory to scrutinize the physics of the strong interactions (for a broad overview see Ref. [1]). In particular, studying their production in relativistic heavy-ion collisions promises vital insight into the properties of quarks and gluons under the extreme conditions of high temperature and density which are present shortly after the Big bang (for a recent review see Ref. [2]). In contrast to the earliest moments of the universe, a heavy-ion collision produces energetic debris that are not in static thermal equilibrium (for selected perspectives see, e.g., Refs. [3][4][5][6][7] and references therein). At the time of incident, hard partonic scatterings may convert the vast amounts of kinetic energy of the incoming projectiles into particles with high velocity or large masses (compared e.g. to the characteristic scales of QCD Λ QCD /p Q 1). It is here where the constituents of quarkonium particles are born. On the other hand, the lighter particles, which are created for example from the fragmentation of the strong initial glasma color fields, quickly thermalize locally and form a liquid-like and expanding quark-gluon plasma. One challenge in this field of study is to understand how the hard particles, which are produced in the initial stages, propagate in the presence of the quasi-thermal hot environment formed by the light degrees of freedom.
Quarkonium in equilibrium with its surroundings has been studied thoroughly in the past using lattice QCD simulations and effective field theories (for some recent works see [8][9][10][11][12]), potential models (see e.g. [13,14]), QCD sum rules [15] and holography (see e.g. [16][17][18]). In recent years, thanks to advances in real-time methods, the focus has shifted into the realm of non-equilibrium physics, heavily motivated by the phenomenologically relevant setting of relativistic heavy-ion collision. A fruitful exchange of ideas between the high-energy nuclear physics community and the condensed matter physics community (see e.g. [19][20][21]) has given momentum to the ongoing development of an open quantum systems description of heavy quarkonium (for the most recent works see e.g. [22][23][24][25][26][27][28][29][30][31][32][33]) in contact but not necessarily in equilibrium with its hot environment (for a recent review see [34]). A central ingredient in the development of these real-time descriptions is the inherent separation of scales between the heavy-quark rest mass and the other characteristic scales, such as energy density and the QCD scale Λ QCD . In vacuum, this separation of scales has helped to establish that the physics of heavy quarkonium bound states may be captured reliably by a non-relativistic potential description. Formalized in the effective field theories Non-relativistic QCD (NRQCD) and potential NRQCD (see [35] for a review), it has been established, at least in a weakly-coupled context, to what extent such a potential picture is also applicable at finite temperature [36][37][38]. The determination of the interaction potential non-perturbatively using lattice QCD simulations is an active field of ongoing research [39][40][41][42][43][44]. Recent work on quarkonium as open-quantum system has shown how the real-and imaginary part of this in general complex valued static heavy quark potential govern the evolution of quarkonium states in a hot medium [27,32]. While the in-medium real part informs us about the screening of the interaction between the heavy quark-antiquark-pair, the imaginary part encodes how scattering with gluons of the surrounding medium over time leads to color decoherence.
The complex static interquark potential at finite temperature thus plays a vital role in the real-time description of quarkonium bound states. It is now understood that we may access its values non-perturbatively, by inspecting the rectangular real-time Wilson loop W (t, r) = Tr Pexp ig which resides on the path that a pair of static color sources traces out as it evolves in real-time. A µ refers to the gauge field of the strong interactions, and, as a correlation function, the Wilson loop is evaluated in path-ordered fashion, indicated by the operator P. Its gauge invariance follows from taking the color trace, denoted above by Tr. As was shown in detail in [45], if the time evolution of a pair of color sources, described by W (t, r), proceeds according to a static Schrödinger equation, we may use it to define the corresponding potential via In the context of the present paper, we consider times bigger than the time corresponding to multiple gluon mediated exchanges as late times. These gluon mediated exchanges may at late times then be collectively replaced by an instantaneous interaction potential. In practice, evaluating eq. (2) often presents technical difficulties.
In turn it has been worked out (see e.g. Ref. [46]) that the values of the potential may be reliably and efficiently extracted from the spectral function of the Wilson loop instead. The spectral function is the unique real-valued function describing the different incarnations of quantum field theoretical correlation functions such as of time ordered, retarded, or Matsubara type. If a time independent potential emerges at late times according to (2), one can show that there must exist a lowest lying peak structure of skewed Breit-Wigner form in the spectral function, whose position and width encode the real and imaginary part of the potential respectively. Around the peak maximum, we may describe it with the following functional form The values of the static potential have been investigated in the past through conventional lattice QCD simulations. However, the fact, that these simulations are carried out in Euclidean time, necessitates the solution of an ill-posed inverse problem to extract the spectral function, often attacked using methods of Bayesian inference (see e.g. [47]). While robust estimates of the real part of the potential have been obtained in this fashion, access to the imaginary part is still severely limited. The reason is that the determination of a spectral width requires significantly higher input data quality than that of spectral peak positions.
In this study, we set out to investigate the complex interquark potential in a genuine real-time setting. Since the notorious sign problem of the quantum path integral prohibits its direct numerical evaluation in Minkowski-time, we will have to agree to a compromise. In case of this study, it amounts to neglecting the quantum fluctuations in the theory and simply focusing on the statistical fluctuations, wherein such statistical fluctuations may be introduced through a thermal medium. I.e., we will resort to the classical statistical approximation of Yang-Mills theory to investigate the binding properties of a pair of static color charges at finite temperature. The classical statistical treatment of gauge fields has a long history in the context of research on Baryogenesis (see e.g. Ref. [48] and references therein), reheating in the early universe in Refs. [49][50][51][52] and more recently in the study of the overoccupied glasma in the early stages of relativistic heavy-ion collisions ( for a review see Ref. [4]). The static interquark potential has been studied in the classical statistical approximation first in Ref. [53]. Here, we will extend and improve on those results, based on work presented in A. Lehmann's doctoral thesis [54].
While the non-equilibrium properties of several phenomenologically relevant systems can be captured quantitatively by classical statistical simulations, we can only expect to gain qualitative insight in the case of thermal equilibrium. The reason is the classical Raleigh Jeans divergence. Thermal modes saturate the spectrum of any classical statistical computation and in turn the effects of, e.g., charge screening, dominated by modes around the cutoff, become dependent on this UV cutoff.
Viewed through the lens of lattice perturbation theory, one furthermore predicts [55][56][57] that the Debye mass changes with the square root of the temperature at fixed lattice spacing which is in contrast to the full quantum theory, where the Debye mass is proportional to the temperature itself. The absence of a well-defined continuum limit precludes us, e.g., from attaching physical units to the simulations. There exist prescriptions of how to make possible a continuum limit by amending the classical simulation with a perturbative treatment of modes close to and above the lattice cutoff (see e.g. [58,59]). For our purposes of qualitatively identifying the form of the interaction potential, the naive implementation suffices, however. Keeping these caveats in mind, classical statistical simulations, as one of the very few direct real-time methods available, promise vital insight into the dynamics of binding of color sources in Yang-Mills theory and provide guidance to understand the fully quantum theory so far only accessible through the lens of static Euclidean lattice QCD simulations.
In the following section 2, we discuss the differences in treatment of static sources between conventional Euclidean lattice QCD simulations and those carried out in the classical statistical approximation. Section 3 briefly describes the numerical methods used in our study before we present our findings in section 4. The paper closes with a conclusion and outlook in section 5.

Static Sources in Lattice Gauge Theory
The study of the binding properties of static charges in the presence of a medium of light charge carriers goes back to the works of Debye and Hückel [60]. From general thermodynamic considerations and use of linear response theory, they deduced classically that in such a system the interactions between the static charges will be screened, i.e., the long-ranged Coulomb interaction will become short ranged with a characteristic screening radius given by the temperature dependent density of light charge carriers. The study of the real-part of the interquark potential in conventional lattice QCD supports this conclusion also for the strong interactions, showing clear signs for a bound state sustaining real-part of the potential in vacuum and a gradual weakening of its values towards a screened form at high temperatures [40].
What would be the intuitive expectation for a classical theory of strong interactions? In the → 0 limit, Yang-Mills theory may be seen as a non-linear extension of Maxwell's electrodynamics. Since charges interact in electrostatics, there is no a priori reason to suspect otherwise in the non-linear theory. Irrespective of the sign of the interaction, we would expect that its strength is encoded in the finite real-part of a classical interaction potential.
The results of Ref. [53] hence are puzzling. The authors set out to study the static potential in classical statistical lattice gauge theory based on the real-time Wilson loop introduced in section 1. They observe that its values are purely real. According to eq. (2), if the Wilson loop goes over into a single exponential behavior at late times it depends on the values of the potential as A purely real W (r, t) thus corresponds to a purely imaginary V (r). Does this mean that the real-part is screened to such small values that it is practically invisible? Or does it mean that there are no interactions between color sources present at all in the classical statistical theory? In this study, we propose and show that the potential between static color sources is complex. It possesses a non-vanishing real and imaginary part. The previously observed absence of a real part originates from subtle differences in how static sources have to be treated in the quantum versus the classical theory.
In fully quantum lattice gauge theory, the physics of static color sources are conventionally investigated by computing the expectation values of the imaginary time Wilson loop W (−iτ, r) in the Euclidean path integral. At first, simulations of the gauge (and possibly also light fermion) degrees of freedom are carried out, which are oblivious of the physics of the static sources that we wish to investigate. It is the evaluation of the Wilson loop itself that introduces these charges into the system. How is this accomplished? To not encumber us with technical difficulties, let us assume that we are in spatial axial gauge, so that only the temporal stretches of the Wilson loop remain relevant: The matrix M denotes a 3 × 3 matrix that arises when absorbing the two exponentials into the exponent of the Feynman weight, forming an effective action in the presence of sources. We see that introducing the Wilson loop as observable actually amounts to a reweighting of the simulated system without static sources to a system with static sources present at exactly those spatial positions, at which the temporal stretches of the Wilson loop start and end. It is this possibility to introduce a posteriori the sources into the system, that makes conventional lattice QCD simulations so versatile. One set of numerical simulations can be reused to investigate multiple different physics scenarios, simply by a choice of observable. The situation is quite different in the classical statistical approximation. There, we do not evaluate the path integral itself but instead evolve fields via deterministic classical equations of motion, while their initial values are drawn from a statistical ensemble. The first difference, we must note, is that computing an observable on the field configurations evolved in that way is not the same as evaluating the observable inside the full path integral. Hence, evaluating the Wilson loop on the classical equations of motion does not amount to the same reweighting we have seen taking place in the full quantum theory. I.e., evaluating the Wilson loop here leaves the gauge fields still oblivious of the presence of the static sources.
As has been shown in detail in [61], the classical equations of motion emerge naturally in the derivation of the classical statistical approximation from the full path integral. We may write the classical limit of the partition function for the gauge fields in the presence of fermions, denoted by Q(x), as Z src in the following form What is left of the full path integral is an integration over the statistical distribution ρ of the classical field A and its conjugate momenta Π at initial time. The dynamics of the fields at later times are determined by the classical equation of motion, housed in the delta function term, featuring the covariant derivative D µ = ∂ µ − igA µ . The coupling of fermionic degrees of freedom to the gauge fields in the classical statistical approximation introduces a current j a µ , coupled to the classical gauge field degrees of freedom. This current is intimately related to the fermion propagator The gauge coupling g appears together with the generators of the gauge group T a and the gamma matrices γ µ , required for coupling the fermions to the gauge field. For heavy fermions, the states, over which the expectation value is taken, do not contain any heavy quark particles. Thus, the commutator reduces to the forward correlation function where only creation operators act on the states from the left. Applying the Foldy-Tani-Wouthuysen transform to separate the contents of the four component Dirac spinor Q = (ψ, χ) into two two-component Pauli spinors, one finds that only the zeroth component of the current remains relevant. Its expression further simplifies to All other spatial components of this current are suppressed and therefore the information of the static sources only enters via the color density j a 0 . It makes its appearance in the classical theory as a source term on the RHS of the Gauss's law written here concisely using the color components of the covariant derivative D µ = ∂ µ − igA µ and the field strength tensor F µν = [D µ , D ν ]. Being interested in the physics of color singlet states, we initialize the system with two static sources placed at x 1 and x 2 and select a color combination encoded in the 3 × 3 matrices M 1 and M 2 . We arrive at the proper Gauss's law for the gauge fields in the presence of two static sources in the classical statistical approximation The matrix' entries may be ordered as "red, green, blue". For a color-anticolor structure relevant to describe an overall color singlet state, we may then choose a red anti-red rr-configuration, which we would represent by M 1 = diag[+1, 0, 0] and M 2 = diag[−1, 0, 0]. This proper Gauss's law, as part of the equations of motion, implements the back-reaction of the sources onto the gauge fields. In turn, missing the source term amounts to neglecting this back-reaction all together. As we will see in the following sections, this fact allows us to understand the absence of a real part in the potential observed in Ref. [53]. In that study, the propagating static color sources, described by the Wilson loop, did not lead to a back-reaction onto the gauge fields, which in turn did not allow them to build up a force among them. On the other hand, scattering of the medium gluons with the heavy quarks does not require such a back coupling and thus a finite imaginary part was found.
In the present study, we will carry out simulations of the Wilson loop in classical statistical gauge theory in the presence of sources, i.e., based on the discretized counterpart of the proper Gauss's law of eq. (11). It allows us to reveal the presence and values of the real part of the static potential. In addition, we investigate whether the back-reaction changes the behavior of the imaginary part of the potential found in Ref. [53].

Numerical Methods
Let us briefly recollect the standard techniques of classical statistical lattice gauge theory for the gluon field degrees of freedom. Our starting point is the anisotropic Wilson plaquette action [62,63] in the presence of a static charge density j a 0 (x). We place it on a hypercubic lattice with spatial and temporal lattice spacing a s and a t respectively The plaquettes are defined as elementary Wilson loops , which contain the discretized gauge field A a µ (x) and the generators of SU (3), denoted by T a . In the classical theory, the gauge coupling g can be absorbed into a redefinition of the gauge fields and will therefore be set to unity in the following. The plaquettes are related to the gauge field strength tensor P µν (x) = exp ia µ a ν F a µν T a ] and in turn allow us to explicitly define the electric fields as In order to simulate this system in a reliable fashion, we go over to a Hamiltonian formulation, where only space remains discretized and time becomes continuous. The strategy aims at formulating the equations of motion in canonical form, which in turn enables the deployment of symplectic time-stepping prescriptions, such as the O(∆t 2 ) accurate leap-frog algorithm for separable Hamiltonians. The canonical form of the equations of motion are obtained by choosing a suitable gauge. In temporal gauge with U 0 (x) = 1, the dynamical degrees of freedom are the gauge group valued spatial links U i (x) and their generator valued conjugate momenta E a i (x), residing on individual time slices. The Hamiltonian according to the plaquette action reads where in the last line we have introduced the dimensionlessH, explicitly recovering that the units of H are that of the inverse spatial lattice spacing. The dynamics of the gauge fields is described via the spatially discretized e.o.m.
Here, S j and S j denote the backward and forward staple, which are the leftovers of the full backward and forward plaquettes where the link U j has been removed, i.e., they yield the backward and forward plaquette once multiplied with the link U j from the right or left respectively. In addition to these equations of motion, the functional derivative of the lattice action with respect to the A 0 component of the gauge field leads to the Gauss's law in the presence of static sources If consistent initial conditions are chosen such that they fulfill eq. (16), then the evolution equations eq. (15) will preserve this property. When implemented as discrete time stepping, the dynamics will lead to small deviations from Gauss's law, which we check to remain insignificant at all times. In order to generate the collection of stochastic initial conditions according to the classical density matrix at temperature T = 1/β, we resort to the strategy originally devised in [64][65][66] and its implementation as outlined in e.g. [67].
The thermal probability distribution of the degrees of freedom follows the Boltzmann weight P [E, U ] ∝ exp − βH[E, U ]]. As the Hamiltonian is separable, the electric field contribution is independent of the links and it appears at first sight that the E a i s are simply Gaussian distributed. Gauss's law however distinguishes among the electric fields those that are physical and those that are not. Hence we will have to not only stochastically draw values of E a i according to but in addition project to the physical subspace of electric fields by subsequently minimizing, via gradient descent, the functional x,a Tr[T a G(x, t)] 2 . Note that it is only here in the projection of the initial conditions where the effect of the static sources enters the gauge field dynamics, as it is here where the Gauss's law is enforced.
With a first set of projected electric fields at hand, we may proceed to compute the corresponding spatial links. To this end, we evolve some arbitrary initial set of links according to eq. (15) based on the quasi-thermal electric fields found above. This allows them to mutually equilibrate with a temperature, which however is not yet the desired one. To reach the thermal steady state between E a i and U i , defined by β, the electric fields need to be redrawn, projected and evolved with the links several times. The outcome of this procedure in turn is used as one of the realizations of the thermal initial conditions for the actual real-time evolution of the gauge fields. Since the steady state is characterized by a constant energy, we have made sure that the mutual equilibration between links and electric fields is realized well enough, so that any remaining changes in the energy are below percent level.
In temporal gauge the computation of the Wilson loop simplifies considerably, as the temporal links are unit matrices. Due to time translational invariance in thermal equilibrium, we compute the real-time Wilson loop spanning from the initial time slice to the time slice, which the simulation has currently reached. To this end we keep a copy of the initial conditions in memory and form the following products In case that the sources are set to zero, the Wilson loop expectation value does not depend on the axis direction and position of it starting point so that W no src (r, t) = x j W j x (r, t) / x j 1. On the other hand in the presence of sources, we choose to place one of them at the origin and the second one along the x-axis at distance r src , such that the Wilson loop at only a single position and with a single spatial extent is computed over time W src (r src , t) = W x 0 (r src , t) . The lack of geometric averaging in this case leads to significantly larger statistical uncertainties in the observable, which in turn need to be compensated for by collecting more statistics for the ensemble average over initial conditions.
As discussed in section 1, the interaction potential between the static charges may be computed using the spectral function of the Wilson loop. While in standard Euclidean lattice simulations this requires to solve an ill-posed inverse problem, in a real-time simulation we only have to carry out a Fourier transform. To minimize the computational cost, we exploit that the Wilson loop at positive and negative times is related by complex conjugation W (t, r) = W * (−t, r) in order to carry out a discrete Fourier transform on its values in an interval [−t max , t max ]. Since the discrete Fourier transform is defined for periodic signals, it is paramount that the simulation has progressed far enough so that the values of the Wilson loop at the latest time have decayed close enough towards zero. This is realized in most cases except for the smallest spatial separation distances which may lead to ringing artefacts in the resulting spectral functions. The further the amplitude of the signal has decayed over time the weaker the ringing becomes. As is well known, such ringing can be taken care of by introducing appropriate windowing functions. We deploy the well established Hann window function [68] to the real-time correlation functions at the lowest and next to lowest spatial separations in order to avoid such ringing artifacts. We have checked that this procedure does not introduce a bias onto the estimation of the potential beyond the statistical uncertainties.
With the spectral function of the Wilson loop at hand we are in a position to investigate the binding properties of static color charges in lattice Yang-Mills theory in classical thermal equilibrium.

Numerical Results
The main results of our study are computed based on lattices with N = 32 grid points in each spatial direction (The source code for our simulations is freely available under an open access license, hosted at the Zenodo repository [69]). We set the spatial lattice spacing to unity a s = 1 and work with a time stepping ∆t = 0.1. Each realization of the microscopic classical dyanmics, started from an independent stochastically drawn set of N runs initial conditions, is evolved up to times N t ∆t = 1200, both in the absence and in the presence of sources. To generate the thermal initial conditions we draw and project the electric fields N init = 20 times and each time evolve them with the gauge fields for N therm = 150 steps using the above ∆t. Supplementary plots showcasing the initialization as well as the conservation of the Gauss's law and energy can be found in appendix A.  Table 1 Parameters used in the simulations underlying our extraction of the proper static interaction potential in classical lattice gauge theory.

Simulations without explicit source term
As a first step, we reproduce and refine the results of [53], in which the thermal realtime Wilson loop was studied in the absence of an explicit source term in Gauss's law. In order to resolve the changes due to variation in the system temperature, we choose ten values of β = 1/T a s between 4 and 32. Using slightly different conventions than in [53], our choice of β = 4 corresponds to β L = 16 in the work of Laine and Tassler. As a representative example, we plot in fig. 1 the Wilson loop's expectation value at β = 4 along real-time t evaluated at several different spatial distances r/a s = 1 . . . 10 (darker to lighter lines). β = 4 corresponds to the highest temperature considered in this study and the falloff of the Wilson loop is thus the most rapid, requiring 500 independent initial conditions to reach the level of statistical uncertainty, indicated by the errorbands.   The top panel of fig. 1 displays the real-part of W no src (r, t) , while the bottom panel shows its imaginary part. As expected from Ref. [53], the Wilson loop in the absence of static sources is purely real within statistical uncertainty. A first inspection of this logarithmic plot by eye indicates that at late times it decreases approximately as a single exponential. It is this exponent that has previously been identified with the imaginary part of the potential between static quark-anti-quark sources.
In order to extract the potential from the Wilson loop, we proceed to scrutinize its spectral function, computed via discrete Fourier transform, plotted as individual datapoints in fig. 2 for six different spatial separation distances r/a s = 2, 4, 6, 8, 10 and 12. We plot as range of frequencies ωa s ∈ [−0.25, 0.25], where the relevant structure encoding the potential is present. The full frequency range extends to the much larger values given by π/∆t (not shown). In agreement with fig. 1, the spectral function is purely real and and exhibits a single dominant peak located at vanishing frequency, whose width increases with spatial separation distance of the underlying coordinate space Wilson loop. The solid lines correspond to the best fit according to eq. (3), using 15 of the discrete points of the spectral function above and below the peak maximum. The skewed Breit-Wigner shape predicted on general grounds in [45] matches the data very well with χ 2 /d.o.f. ≈ 1.5 . . . 2. Changes in the fitting range do not significantly change its outcome. In order to assign uncertainties to the extracted values of the potential we have deployed a ten bin Jackknife, whose variance underlies the relatively small errorbars shown in fig. 3.
In order to crosscheck the results obtained from the spectral function analysis, we have also carried out exponential fits to the late time behavior of the coordinate space Wilson loop. We see that at small separation distances r/a < 6 where the single exponential decay of the Wilson loop is well resolved, both methods agree. At larger spatial separation distances, it becomes more and more difficult to select by eye an appropriate fitting regime, leading to exponential fit values for ImV that are slightly larger than those from the spectral function analysis, an effect we attribute to excited state contamination. Considering the agreement with the exponential fit at small distances and the fact that the DFT based extraction shows very small dependence on the fitting range, we deem it reliable enough in the context of the currently available input data quality.
Tabulating the values of the peak position and width for different separation distances provides us with an estimate of the static quark interaction potential in the absence of a back-reaction of those sources onto the gauge fields. Carrying out the same procedure at different values of β, we elucidate the temperature dependence of the potential which is plotted in the top panel of fig. 3. What we obtain from the width are positive numbers. These correspond to the magnitude of the imaginary part. Its values are surely negative, as they lead to a dampening of the amplitude of the Wilson loop over time. With the imaginary part in the full quantum theory possessing a trivial dependence on temperature, we plots here its values divided by the system temperature. We find an imaginary part, which increases monotonously with temperature. Based on the state-of-the-art extraction procedure using the Wilson loop spectral functions, the values of Im[V ] obtained here are robust at the current level of statistical uncertainty up to r/a s = 16 for β ≥ 16, while the loss of signal to noise ratio limits the extraction to smaller distances for lower values of β.
Qualitatively similar to the predictions of HTL perturbation theory in the fully quantum case, we find two distinct regions. One, at small distances (here r/a s < 8), which features a relatively steep increase and a convex shape. The other at subsequently larger distances, with a much weaker slope and a concave behavior that appears to lead to a flattening of the values of ImV at larger distances.
We may compare the simulation results also quantitatively to the predictions from the classical limit of lattice perturbation theory. Leveraging equation (3.23) from Ref. [53], we compute the values of |ImV |/g 2 T for the β values in our study (see the supplementary material for an explicit implementation in Mathematica). The outcome is plotted in the top panel of fig. 3 as gray solid lines. Laine and Tassler evaluated the corresponding β = 4 imaginary part for the first four spatial distances. We see that while the qualitative behavior appears similar, quantitative differences exist. In general the perturbative values are smaller than those obtained from the lattice. In addition the perturbative ImV reaches its asymptotic value at earlier distances compared to what we observe in the simulation. It is not surprising that also the large distance asymptotic value is smaller in the weakly coupled case, as ImV(r → ∞) is related to single heavy quark energy loss in the medium. A strongly coupled environment intuitively induces energy loss more efficiently than a weakly coupled one.

Simulation with explicit source terms
In this section we advance toward the main result of our study, the determination of the proper interaction potential between static color sources in classical lattice gauge theory at finite temperature. As discussed in section 2 the missing ingredient in studies so far was the inclusion of the back-reaction of the static sources on the gauge fields. Now with source terms explicitly present, we have to carry out separate simulations for each of the different spatial separation distances we wish to observe the static color sources at. In addition, it is not possible to average the Wilson loop over axis-directions and spatial positions, which requires a factor of up to 10 times higher statistics in order to reach satisfactory results. We therefore restrict the study to one single temperature given by β = 20. Our choice is motivated by the fact that at this β, in simulations without explicit source terms present, thermal effects are already visible, while at the same time a sufficient signal-to-noise ratio can be obtained.
As discussed in section 3, the novel element in this study is the treatment of the thermal initial conditions in order to take into account the change in Gauss's law due to the sources. We have checked explicitly (see appendix A) that this proper Gauss's law is preserved by the dynamical evolution.
In fig. 4, we show the resulting real-(top) and imaginary part (bottom) of the Wilson loop at different spatial separation distances r/a s = 1 . . . 10 (darker to lighter solid lines) at inverse temperature β = 20. The expectation values are computed from N runs = 1200 simulations, carried out with independently drawn initial conditions. The simple modification in the initial conditions has fundamentally changed the behavior we observe.
Not only does the imaginary part of the Wilson loop now show finite values, Re[V ] and Im[V ] both exhibit oscillatory behavior with a monotonously decreasing amplitude also. Note that the former starts from unity while the latter starts from zero at the origin, consistent with the transformation properties of the Wilson loop under time reversal. A first inspection by eye hints at the presence of at least two sinusoidal contributions to the dynamics. The exponential suppression of the amplitude of the Wilson loop with increasing spatial separation distances is also present here, which leads to a quickly deteriorating signal to noise ratio at increasing r/a s .   While the decrease in amplitude hints at the presence of a finite imaginary part of the potential, the oscillatory component bodes well to identify a real-part of the interaction potential as well.
Let us proceed to a quantitative extraction of the potential for which we first compute the spectral function of the Wilson loop via DFT. In fig. 5 we plot the corresponding values (colored symbols) for positive frequencies in the limited domain between ωa s ∈ [0.045, 0.12] where we locate one of the two dominant features of the spectrum, which is directly related to the potential. The results for six spatial separation distances are plotted between r/a s = 1 and 6 (darker to lighter color).
Due to the significantly lower signal to noise ratio, the data points here show much more variation than in the case without sources. Since the amplitude of the r/a s = 1, 2 Wilson loop has not yet decayed significantly at the final time ta s = 1200 of our simulation, the naive Fourier transform would exhibit ringing around the lowest lying peak. We thus have applied the Hann windowing function in time to the coordinate space Wilson loop before the DFT. We find that also in the case with sources present, the lowest lying spectral structure, around its maximum, can be captured well with the skewed Breit-Wigner form of eq. (3), as shown by the best fit results as solid lines. Both the peak position and the peak width show a clear dependence on the separation distance. By determining the peak positions and widths at β = 20, we arrive at the central result of this study, the values of the real-and imaginary part of the proper static interquark potential, shown in Figure 6. As expected by the presence of a finite imaginary part in the coordinate space Wilson loop W src (r, t) and its damped oscillatory behavior, we do indeed find finite values for the real-part, as shown in the top panel of fig. 6. The errors are dominated by the fit uncertainty of the sevenparameter χ 2 fit according to eq. (3), which has been carried out for ten Jackknife bins.
Since classical lattice gauge theory does not capture the physics of confinement, one may expect to find a Coulombic behavior at very small distances corresponding to the tree-level short range interaction in Yang-Mills theory. Since we are observing the system at finite temperature, a hot medium of color charges fills the space between the static sources and in analogy with the well-known Abelian theory, we thus expect that the interactions are screened. And indeed using the simple ansatz − A r exp[−m D r] + const. it is possible to reproduce the dependence of Re[V ] on the spatial separation distance r very well (dashed green line). Such a Debye screened behavior is qualitatively very similar to what has been observed also in the fully quantum theory deep in the deconfined phase.
The value of m D we observe is in good agreement with that predicted by lattice perturbation theory according to eq.  In the bottom panel of fig. 6, we plot the absolute value of the imaginary-part of the proper potential as colored symbols, again divided by the temperature. A very similar picture emerges as in the simulations without explicit sources. The imaginary part increases monotonously with separation distance. Plotting the results from simulations without explicit inclusion of sources at the same temperature as light gray symbols, we actually find very good agreement at those distances where the extraction in the presence of explicit back-reaction is reliable. We thus reconfirm the results of [53] for the imaginary part, observing no significant effect of the backreaction on the physics of scattering between static sources and medium constituents in the value of Im [V ].
Interestingly, the values of the imaginary part, we obtain in lattice simulations, are significantly larger than those predicted by lattice perturbation theory (compare to fig. 3). On the other hand, the Debye mass, governing the behavior of the real-part comes out at a very similar value. It would be interesting to explore, whether the agreement for Im [V ] improves significantly when extending the perturbative calculation of Ref. [53] to one higher order.

Conclusion
Classical statistical simulations of gauge fields in the presence of static sources provide us with qualitative insight into the dynamical binding mechanisms of heavy quarks under the strong interactions. In this paper, we have revisited how the presence of static sources affects the classical equations of motion, leading to an explicit source term in Gauss's Law, which was absent in previous works.
In order to scrutinize the binding of color sources, we extract the static potential acting between them, based on the state-of-the-art approach of computing the Wilson loop spectral function and carrying out a skewed Breit-Wigner fit to its lowest lying structure. In the absence of back-reaction this method allows us to confirm the results of a finite imaginary part first presented in [53] and extend them to larger spatial separation distances and a broader temperature range.
By employing the proper Gauss's Law in the presence of static sources, we observe significant changes in the dynamics of the Wilson loop already on a qualitative level. Its values become complex and show oscillatory behavior which directly translate into finite values of the real-part of the static potential. We investigate its behavior at β = 20 and find that it is well described by a Debye screened form. The imaginary part of the potential interestingly appears not to be affected by the inclusion of the back-reaction when compared to a simulation where the explicit back-reaction through the Gauss's law was absent.
The Wilson loop in conventional gauge theory embodies a reweighting from a theory without sources to the theory of interest with sources present. In the classical statistical theory it does not play this role and one may think about whether it is the most appropriate quantity to extract the interaction potential there. One argument to consider is that in contrast to SU (3) gauge theory deployed in this study, we could just as well take SU (2). In the fully quantum theory both theories are known to lead to interactions between color sources. However since the Wilson loop in SU (2) must be purely real, it being the trace over the product of link variables, it cannot encode a finite real part of the potential, in the form observed here.
An alternative approach to determining the interaction potential has been put forward in the context of Euclidean lattice gauge theory in Ref. [70]. The authors propose to compute the expectation values of the gauge invariant energy momentum tensor and to inspect its spatial components. This stress tensor encodes the force acting locally on the faces of a given volume element. Integrating the net force up, allows one to compute the underlying potential. While this approach has been implemented with success in the Euclidean theory, the fact that the current discretizations of the energy momentum tensor are not conserving in Minkowski time prevent us from utilizing it straight away. It actually requires the development of a new conserving discretization scheme, which is currently work in progress.
Having elucidated the physics of static sources in the presence of classical statistical gauge fields in this study, we set out to compute the interaction between very heavy but still dynamical quarks in classical statistical lattice gauge theory in the future using a real-time implementation of the effective field theory NRQCD (for preliminary work in this direction see [71]  that thermal equilibrium between electric fields and links was reached and that the leap-frog time-stepping as symplectic solver realizes its potential to preserve energy in the discretized e.o.m..
The last check to be made in classical statistical simulations of Yang-Mills fields is whether the discretized dynamics preserve the Gauss's law. And while the time stepping introduces a finite deviation as shown in the bottom panel of fig. 8, it remains fully insignificant over the whole time period of the simulation. Note that the starting value of around 5×10 −13 arises from our choice of tolerance in projecting the electric fields to their physical values.
The simulations in the presence of explicit sources show similarly robust behavior. In fig. 9 we show the change in total energy over the N init = 20 thermalization cycles after which the asymptotic value is well established.
Taking a look at the evolution of the relative change in total energy during the actual simulation in the top panel of fig. 10, we find again that any residual change is on a below the per-mille level and thus as insignificant as in the case without sources. The proper Gauss's law, the vital new ingredient in the determination of the proper interaction potential between static sources in lattice gauge theory, also turns out to be well preserved over the whole simulation time. As shown in the bottom panel of fig. 10 it remains below 10 −13 at all evolution steps. Again the start value of around 4 × 10 −14 is related to the tolerance of projecting the initial electric fields to their physical subspace.