# Proper static potential in classical lattice gauge theory at finite T

ALAlexander Lehmann  ARAlexander Rothkopf  Institute for Theoretical Physics,
Heidelberg University, Philosophenweg 16, 69120 Heidelberg, Germany  Faculty
of Science and Technology, University of Stavanger, 4021 Stavanger, Norway

###### Abstract

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.

static quark potential,

classical lattice gauge theory,

quarkonium,

###### keywords:

Research Article - PreprintFP

##  1 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. Brambilla:2010cs ). 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. Rothkopf:2019ipj ). 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. Strickland:2020hyg ; Berges:2020fwq ; Shen:2020gef ; Akamatsu:2020lej ;
Kurkela:2019kdu  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 $\Lambda_{\rm QCD}/p_{Q}\ll
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 Larsen:2019zqv ; Offler:2019eij ; Larsen:2019bwy ;
Kim:2018yhk ; Burnier:2017bod ), potential models (see e.g. Burnier:2015tda ;
Burnier:2016kqm ), QCD sum rules Song:2020kka  and holography (see e.g.
Braga:2017bml ; Grigoryan:2010pj ; Fujita:2009ca ). 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. Borghini:2011ms ; Akamatsu:2011se ;
Akamatsu:2012vt ) has given momentum to the ongoing development of an open
quantum systems description of heavy quarkonium (for the most recent works see
e.g. Brambilla:2020qwo ; deJong:2020tvx ; Yao:2020eqy ; Yao:2020xzw ;
Alund:2020ctu ; Miura:2019ssi ; Brambilla:2019tpt ; Yao:2018sgn ;
Blaizot:2018oev ; Brambilla:2017zei ; Kajimoto:2017rel ; Yao:2017fuc ) in
contact but not necessarily in equilibrium with its hot environment (for a
recent review see Akamatsu:2020ypb ). 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 $\Lambda_{\rm 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 Brambilla:2004jw  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 Laine:2006ns ; Beraudo:2007ky ;
Brambilla:2008cx . The determination of the interaction potential non-
perturbatively using lattice QCD simulations is an active field of ongoing
research Rothkopf:2011db ; Burnier:2014ssa ; Burnier:2016mxc ;
Petreczky:2017aiz ; Petreczky:2018xuh ; Lafferty:2019jpr . 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 Kajimoto:2017rel ;
Miura:2019ssi . 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

$\displaystyle W_{\square}(t,r)=\left\langle{\rm Tr}\Big{[}{\cal P}{\rm
exp}\Big{[}ig\int_{{\cal
C}_{\square(r,t)}}dx^{\mu}A^{\mu}(x)\Big{]}\Big{]}\right\rangle\,,$ (1)

which resides on the path that a pair of static color sources traces out as it
evolves in real-time. $A_{\mu}$ 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 ${\cal P}$. Its gauge
invariance follows from taking the color trace, denoted above by ${\rm Tr}$.
As was shown in detail in Burnier:2012az , if the time evolution of a pair of
color sources, described by $W_{\square}(t,r)$, proceeds according to a static
Schrödinger equation, we may use it to define the corresponding potential via

$\displaystyle
V(r)=\lim_{t\to\infty}i\partial_{t}W_{\square}(t,r)/W_{\square}(t,r).$ (2)

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.
Burnier:2013fca ) 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

$\displaystyle\begin{split}\rho_{\square}(r,\omega)=&\frac{e^{{\rm
Im}[\sigma_{\infty}](r)}}{\pi}\cdot\\\ &\frac{|{\rm Im}[V](r)|{\rm cos}[{\rm
Re}[\sigma_{\infty}](r)]-({\rm Re}[V](r)-\omega){\rm sin}[{\rm
Re}[\sigma_{\infty}](r)]}{{\rm Im}[V](r)^{2}+({\rm Re}[V](r)-\omega)^{2}}\\\
&+c_{0}(r)+c_{1}(r)t_{Q\bar{Q}}({\rm Re}[V](r)-\omega)\\\
&+c_{2}(r)t_{Q\bar{Q}}^{2}({\rm Re}[V](r)-\omega)^{2}+\cdots.\end{split}$ (3)

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. Rothkopf:2019dzu ). 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. DOnofrio:2014rug  and
references therein), reheating in the early universe in Refs. Micha:2004bv ;
Berges:2008wm ; Berges:2016nru ; Chatrchyan:2020syc  and more recently in the
study of the overoccupied glasma in the early stages of relativistic heavy-ion
collisions ( for a review see Ref. Berges:2020fwq ). The static interquark
potential has been studied in the classical statistical approximation first in
Ref. Laine:2007qy . Here, we will extend and improve on those results, based
on work presented in A. Lehmann’s doctoral thesis phdthesis_AlexanderLehmann .

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 Arnold_1997 ; WATSON_1939 ;
Glasser_2000  that the Debye mass changes with the square root of the
temperature at fixed lattice spacing

$\displaystyle(m_{D}^{\rm lPT})^{2}=2g^{2}TN_{c}\frac{\Sigma_{\rm lat}}{4\pi
a_{s}},\quad\Sigma_{\rm
lat}=\Gamma^{2}(\frac{1}{24})\Gamma^{2}(\frac{11}{24})\frac{\sqrt{3}-1}{48\pi^{2}},$
(4)

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. Bodeker:1998hm ; Bodeker:1999gx ). 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.

##  2 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
huckel1923theory . 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 Burnier:2014ssa .

What would be the intuitive expectation for a classical theory of strong
interactions? In the $\hbar\to 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. Laine:2007qy  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

$\displaystyle\lim_{t\to\infty}W_{\square}(r,t)\sim{\rm
exp}\big{[}-itV(r)\big{]}$ (5)

A purely real $W_{\square}(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_{\square}(-i\tau,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:

$\displaystyle W_{\square}(\tau,r=|{\bf x}_{2}-{\bf x}_{2}|)$ (6)
$\displaystyle=\frac{1}{Z_{\rm no\,src}}\int{\cal D}[A]{\rm Tr}\Big{[}{\rm
exp}\Big{[}ig\int dtA_{0}({\bf x}_{1},t)\Big{]}{\rm exp}\Big{[}-ig\int
dtA_{0}({\bf x}_{2},t)\Big{]}\Big{]}e^{-S_{\rm gluon}[A]}$
$\displaystyle=\frac{1}{Z_{\rm no\,src}}\int{\cal D}[A]{\rm exp}\Big{[}-S_{\rm
gluon}[A]-g\int d^{4}x\,{\rm ReTr}\big{[}A_{0}(x)M\big{(}\delta^{(3)}({\bf
x}-{\bf x}_{1})-\delta^{(3)}({\bf x}-{\bf x}_{2})\big{)}\big{]}\Big{]}$
$\displaystyle=\frac{1}{Z_{\rm no\,src}}\int{\cal D}[A]{\rm exp}\Big{[}-S_{\rm
gluon}[A]-g\int d^{4}x\,{\rm ReTr}\big{[}A_{0}(x)j_{0}(x,{\bf x}_{1},{\bf
x}_{2})\big{]}\Big{]}$ $\displaystyle=\frac{Z_{\rm src}(\tau,r)}{Z_{\rm
no\,src}}.$

The matrix $M$ denotes a $3\times 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 kasper_fermion_2014 , 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_{\rm src}$ in the following form

$\displaystyle Z_{\rm src}$ $\displaystyle=\int{\cal D}A(t=0)\int{\cal
D}\Pi(t=0)\rho(A,\Pi,t=0)\delta\Big{(}D_{\mu}F^{\mu\nu}[A]-j^{\nu}\Big{)}\,.$
(7)

What is left of the full path integral is an integration over the statistical
distribution $\rho$ of the classical field $A$ and its conjugate momenta $\Pi$
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_{\mu}=\partial_{\mu}-igA_{\mu}$.

The coupling of fermionic degrees of freedom to the gauge fields in the
classical statistical approximation introduces a current $j_{\mu}^{a}$,
coupled to the classical gauge field degrees of freedom. This current is
intimately related to the fermion propagator

$\displaystyle j_{\mu}^{a}(x)=\frac{g}{2}{\rm
Tr}\big{[}\langle[\bar{Q}(x),Q(x)]\rangle_{A}\gamma_{\mu}T^{a}\big{]}\overset{T/m_{Q}\ll
1}{=}\frac{g}{2}{\rm
Tr}\big{[}\langle\bar{Q}(x)Q(x)\rangle_{A}\gamma_{\mu}T^{a}\big{]}.$ (8)

The gauge coupling $g$ appears together with the generators of the gauge group
$T^{a}$ and the gamma matrices $\gamma_{\mu}$, 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=(\psi,\chi)$ into two two-component Pauli spinors, one finds that only the
zeroth component of the current remains relevant. Its expression further
simplifies to

$\displaystyle j_{0}^{a}(x)=\frac{g}{2}{\rm
Tr}\big{[}\langle\psi^{\dagger}(x)\psi(x)\rangle_{A}T^{a}-\langle\chi^{\dagger}(x)\chi(x)\rangle_{A}T^{a}\big{]}+{\cal
O}(m_{Q}^{-1}).$ (9)

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

$\displaystyle D^{a}_{i}F_{0i}^{ab}=j_{0}^{b}(x),$ (10)

written here concisely using the color components of the covariant derivative
$D_{\mu}=\partial_{\mu}-igA_{\mu}$ and the field strength tensor
$F_{\mu\nu}=[D_{\mu},D_{\nu}]$. Being interested in the physics of color
singlet states, we initialize the system with two static sources placed at
${\bf x}_{1}$ and ${\bf x}_{2}$ and select a color combination encoded in the
$3\times 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

$\displaystyle D^{a}_{i}F_{0i}^{ab}\overset{m\to\infty}{=}\frac{g}{2}{\rm
Tr}\Big{\\{}T^{b}\big{(}M_{1}\delta^{(3)}({\bf x}-{\bf
x}_{1})+M_{2}\delta^{(3)}({\bf x}-{\bf x}_{2})\big{)}\Big{\\}}\,.$ (11)

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 $r\bar{r}$-configuration, which we would
represent by $M_{1}={\rm diag}[+1,0,0]$ and $M_{2}={\rm 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. Laine:2007qy . 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. Laine:2007qy .

##  3 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 Wilson_1974 ; Klassen_1998
in the presence of a static charge density $j_{0}^{a}(x)$. We place it on a
hypercubic lattice with spatial and temporal lattice spacing $a_{s}$ and
$a_{t}$ respectively

$\displaystyle S[U]=\frac{1}{g^{2}}\sum_{t}a_{t}\sum_{x}a_{s}^{3}\Big{(}$
$\displaystyle\sum_{i}\frac{2}{(a_{s}a_{t})^{2}}{\rm
ReTr}\big{[}1-P_{0i}(x)\big{]}$ $\displaystyle-$
$\displaystyle\sum_{i,j}\frac{1}{a_{s}^{4}}{\rm
ReTr}\big{[}1-P_{ij}(x)\big{]}+A_{0}^{a}(x)j_{0}^{a}(x)\Big{)}$ (12)

The plaquettes are defined as elementary Wilson loops
$P_{\mu\nu}(x)=U_{\mu}(x)U_{\nu}(x+\hat{\mu})U^{\dagger}_{\mu}(x+\hat{\nu})U^{\dagger}_{\nu}(x)$
and are composed of the gauge links $U_{\mu}(x)={\rm
exp}[ia_{\mu}gA^{a}_{\mu}(x)T^{a}]$, which contain the discretized gauge field
$A^{a}_{\mu}(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_{\mu\nu}(x)={\rm
exp}\big{[}ia_{\mu}a_{\nu}F_{\mu\nu}^{a}T^{a}]$ and in turn allow us to
explicitly define the electric fields as
$E_{i}(x)=E^{a}_{i}(x)T^{a}=F_{0i}^{a}(x)T^{a}$.

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 ${\cal O}(\Delta 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

$\displaystyle H$ $\displaystyle=\Big{[}L_{E}-L_{M}\Big{]}$
$\displaystyle=\sum_{x}\Big{\\{}\sum_{i,a}\frac{a_{s}^{3}}{2}\big{(}E^{a}_{i}({\bf
x},t)\big{)}^{2}+\sum_{i,j}\frac{1}{a_{s}}{\rm ReTr}\big{[}1-P_{ij}({\bf
x},t)\big{]}\Big{\\}}$
$\displaystyle=\frac{1}{a_{s}}\sum_{x}\Big{\\{}\sum_{i,a}\frac{1}{2}\big{(}a_{s}^{2}E^{a}_{i}({\bf
x},t)\big{)}^{2}+\sum_{i,j}{\rm ReTr}\big{[}1-P_{ij}({\bf
x},t)\big{]}\Big{\\}}=\frac{1}{a_{s}}\bar{H},$ (13)

where in the last line we have introduced the dimensionless $\bar{H}$,
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.

$\displaystyle\dot{U}_{i}({\bf x},t)=ia_{s}E^{a}_{i}({\bf x},t)T^{a}U_{i}({\bf
x},t),$ (14) $\displaystyle\dot{E}_{i}^{a}({\bf x},t)=\frac{2}{a_{s}^{3}}{\rm
ImTr}\Big{[}\sum_{j\neq i}T^{a}U_{j}({\bf
x},t).(S^{\sqsubset}_{j}+S^{\sqsupset}_{j})\Big{]}$ (15)

Here, $S^{\sqsubset}_{j}$ and $S^{\sqsupset}_{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

$\displaystyle G({\bf x},t)=\frac{1}{a_{s}a_{t}}\sum_{i}\Big{[}E_{i}({\bf
x},t)-U_{i}^{\dagger}({\bf x},t)E_{i}({\bf x},t)U_{i}({\bf
x},t)\Big{]}-\frac{1}{a_{t}}j^{a}_{0}({\bf x})T^{a}=0.$ (16)

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/\beta$, we resort to the
strategy originally devised in Ambjorn:1990pu ; Ambjorn:1997jz ;
Grigoriev:1988bd  and its implementation as outlined in e.g. Akamatsu:2015kau
.

The thermal probability distribution of the degrees of freedom follows the
Boltzmann weight $P[E,U]\propto{\rm exp}\Big{[}-\beta 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_{i}^{a}$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_{i}^{a}$ according to

$\displaystyle
E^{a}_{i}=\eta^{a}_{i},\quad\langle\eta^{a}_{i}\rangle=0,\quad\langle\eta^{a}_{i}\eta^{b}_{k}\rangle=\sigma^{2}\delta_{ab}\delta_{ik},\quad\sigma^{2}=\frac{1}{a_{s}^{3}\beta},\,$
(17)

but in addition project to the physical subspace of electric fields by
subsequently minimizing, via gradient descent, the functional $\sum_{x,a}{\rm
Tr}[T^{a}G({\bf 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_{i}^{a}$ and $U_{i}$, defined by $\beta$, 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

$\displaystyle W^{j}_{\bf x}(r,t)={\rm
Tr}\Big{[}\prod_{l=0}^{r/a_{s}}U_{j}({\bf
x}+a_{s}l\hat{j},0)\Big{(}\prod_{l=0}^{r/a_{s}}U_{j}({\bf
x}+a_{s}l\hat{j},t)\Big{)}^{\dagger}\Big{]}\,.$ (18)

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 $\langle W_{\square}^{\rm no\,src}(r,t)\rangle=\langle\sum_{\bf
x}\sum_{j}W^{j}_{\bf x}(r,t)\rangle/\sum_{\bf x}\sum_{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_{\rm src}$, such that the
Wilson loop at only a single position and with a single spatial extent is
computed over time $\langle W_{\square}^{\rm src}(r_{\rm
src},t)\rangle=\langle W^{x}_{\bf 0}(r_{\rm src},t)\rangle$. 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_{\square}(t,r)=W^{*}_{\square}(-t,r)$ in order to carry out a discrete
Fourier transform on its values in an interval $[-t_{\rm max},t_{\rm 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 Hann_window  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.

##  4 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
alexander_rothkopf_2020_4332406 ). We set the spatial lattice spacing to unity
$a_{s}=1$ and work with a time stepping $\Delta t=0.1$. Each realization of
the microscopic classical dyanmics, started from an independent stochastically
drawn set of $N_{\rm runs}$ initial conditions, is evolved up to times
$N_{t}\Delta 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_{\rm init}=20$ times and each time evolve them with the gauge fields
for $N_{\rm therm}=150$ steps using the above $\Delta t$. Supplementary plots
showcasing the initialization as well as the conservation of the Gauss’s law
and energy can be found in appendix A.

$\beta$ | $a_{s}$ | src | $N_{\rm runs}$ |
---|---|---|---|---
4 | 1 | NO | 500 |
6 | 1 | NO | 200 |
8 | 1 | NO | 80 |
10 | 1 | NO | 80 |
12 | 1 | NO | 80 |
16 | 1 | NO | 80 |
20 | 1 | NO | 80 |
20 | 1 | YES | 1200 |
24 | 1 | NO | 80 |
28 | 1 | NO | 80 |
32 | 1 | NO | 80 |
Table 1: Parameters used in the simulations underlying our extraction of the
proper static interaction potential in classical lattice gauge theory.

Figure 1: Real- (top) and imaginary (bottom) part of the real-time Wilson loop
evaluated at different spatial on-axis separation distances between $r=1\ldots
10a_{s}$ on $N_{x}=32$ lattices with $\beta=4$ in the absence of a back-
reaction between static sources and the surrounding medium. The real part of
the Wilson loop is positive and exhibits clear exponential damping, while the
imaginary part is compatible with zero.

###  4.1 Simulations without explicit source term

As a first step, we reproduce and refine the results of Laine:2007qy , in
which the thermal real-time 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 $\beta=1/Ta_{s}$
between $4$ and $32$. Using slightly different conventions than in
Laine:2007qy , our choice of $\beta=4$ corresponds to $\beta_{\rm 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 $\beta=4$ along real-time $t$ evaluated
at several different spatial distances $r/a_{s}=1\ldots 10$ (darker to lighter
lines). $\beta=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 $\langle W^{\rm
no\,src}_{\square}(r,t)\rangle$, while the bottom panel shows its imaginary
part. As expected from Ref.Laine:2007qy , 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 $\omega a_{s}\in[-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 $\pi/\Delta 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.

Figure 2: Discrete Fourier transform (colored points) of the real-time Wilson
loop evaluated at different spatial on-axis separation distances between
$r=1\ldots 12a_{s}$ on $N_{x}=32$ lattices with $\beta=4$ in the absence of a
back-reaction between static sources and the surrounding medium. Solid lines
indicate the best fit result of the dominant peak structure at zero
frequencies.

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 Burnier:2012az
matches the data very well with $\chi^{2}/{\rm d.o.f.}\approx 1.5\ldots 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 ${\rm 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.

Figure 3: Real- (top) and imaginary part of the static quark potential
evaluated at different temperatures $\beta=4\ldots 32$ on $N_{x}=32$ lattices
in the absence of a back-reaction between static sources and the surrounding
medium. (bottom) comparison of the results at the four highest temperatures to
the corresponding values from lattice perturbation theory.

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 $\beta$, 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 ${\rm Im}[V]$ obtained here are robust at
the current level of statistical uncertainty up to $r/a_{s}=16$ for $\beta\geq
16$, while the loss of signal to noise ratio limits the extraction to smaller
distances for lower values of $\beta$.

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
${\rm Im}V$ 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. Laine:2007qy , we compute the values of $|{\rm Im}V|/g^{2}T$
for the $\beta$ 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 $\beta=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 ${\rm 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 $\rm{Im}V(r\to\infty)$ 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.

###  4.2 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 $\beta=20$. Our choice
is motivated by the fact that at this $\beta$, 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\ldots 10$
(darker to lighter solid lines) at inverse temperature $\beta=20$. The
expectation values are computed from $N_{\rm runs}=1200$ simulations, carried
out with independently drawn initial conditions. The simple modification in
the initial conditions has fundamentally changed the behavior we observe.

Figure 4: Real- (top) and imaginary (bottom) part of the real-time Wilson loop
evaluated with sources placed at different spatial on-axis separation
distances between $r=1\ldots 10a_{s}$ on $N_{x}=32$ lattices with $\beta=20$.
In the presence of back-reaction, both quantities show finite values and
exhibit a combination of oscillatory and damping behavior.

Not only does the imaginary part of the Wilson loop now show finite values,
${\rm Re}[V]$ and ${\rm 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 $\omega a_{s}\in[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).

Figure 5: Discrete Fourier transform (colored points) of the real-time Wilson
loop evaluated in the presence of sources placed at different spatial on-axis
separation distances between $r=2\ldots 6a_{s}$ on $N_{x}=32$ lattices with
$\beta=20$. For clarity only the positive frequencies of the mirror symmetric
function are shown. Solid lines indicate the best fit result of the lowest
lying peak structure around its maximum, relevant for the extraction of the
static quark potential.

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.

Figure 6: Real- (top) and imaginary (bottom) part of the static quark
potential evaluated in the presence of sources at $\beta=20$ on $N_{x}=32$
lattices. As expected from the behavior of the real-time Wilson loop we now
find both a finite real- and imaginary part. The former exhibits clear signs
of screening as indicated by the Debye fit (dashed gray line). The latter
shows the expected monotonous increase with spatial distance. Results in the
absence of sources are shown as gray symbols and reveal a very good agreement
of ${\rm Im}[V]$ between the simulations with and without back reaction.

By determining the peak positions and widths at $\beta=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
$\langle W_{\square}^{\rm src}\rangle(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
seven-parameter $\chi^{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
$-\frac{A}{r}{\rm exp}[-m_{D}r]+{\rm const.}$ it is possible to reproduce the
dependence of ${\rm 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. 4:

$\displaystyle(m_{D}a_{s})(\beta=20)=0.25(5)$ $\displaystyle(m_{D}^{\rm
lPT}a_{s})(\beta=20)=0.275$ (19)

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
Laine:2007qy  for the imaginary part, observing no significant effect of the
back-reaction on the physics of scattering between static sources and medium
constituents in the value of ${\rm 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 ${\rm Im}[V]$
improves significantly when extending the perturbative calculation of Ref.
Laine:2007qy  to one higher order.

##  5 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 Laine:2007qy  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 $\beta=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.yanagihara_distribution_2019 . 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 Lehmann:2020kjg
).

## Acknowledgements

The authors thank Y. Akamatsu for valuable and insightful discussions. A.L.
was supported in full by the German Research Foundation (DFG) funded
Collaborative Research Centre “SFB 1225 (ISOQUANT)”. A.R. gladly acknowledges
support by the Research Council of Norway under the FRIPRO Young Research
Talent grant 286883. The numerical simulations have been carried out on
computing resources provided by UNINETT Sigma2 - the National Infrastructure
for High Performance Computing and Data Storage in Norway under project
NN9578K-QCDrtX ”Real-time dynamics of nuclear matter under extreme conditions”

##  Appendix A Supplemental figures

In order for the reader to judge the fidelity of our numerical simulations, we
present here representative plots concerning the thermalization procedure, the
conservation of the Gauss’s law as well as conservation of energy. The first
three figures are based on $\beta=4$ simulations in the absence of explicit
sources from the Gauss’s law. The remaining ones refer to the $\beta=20$
simulations in the presence of sources.

Constructing the thermal initial conditions constitutes a major numerical cost
in the classical statistical simulation of gauge fields. As shown in fig. 7
the total energy of the electric fields and spatial links after multiple
mutual equilibration cycles approaches an asymptotic value, characteristic for
that ensemble of thermal initial conditions. In order to keep the cost as low
as possible we have chosen to draw and project the initial electric fields
$N_{\rm init}=20$ times, after which the asymptotic state appears well
realized.

Figure 7: Thermalization of the initial conditions in the absence of sources
at $\beta=4$.

And indeed, when we observe the relative change in the total energy over the
whole dynamical simulation, as plotted in the top panel of fig. 8, we find
that its value remains constant down to the sub-permille level. This in turn
both tells us 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..

Figure 8: Energy conservation (top) and Gauss’s Law preservation (bottom) in
the absence of sources at $\beta=4$.

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\times 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_{\rm
init}=20$ thermalization cycles after which the asymptotic value is well
established.

Figure 9: Thermalization of the initial conditions in the presence of sources
at spatial distance r=a at $\beta=20$.

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.

Figure 10: Energy conservation (top) and Gauss’s Law preservation (bottom) in
the presence of sources at spatial distance r=a at $\beta=20$.

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\times 10^{-14}$ is related to the
tolerance of projecting the initial electric fields to their physical
subspace.

## Competing interests

The authors declare that they have no competing interests.

## Author’s contributions

  * •

A. Lehmann: implementation real-time simulation, partial data analysis,
writing

  * •

A. Rothkopf: conceptualization, funding acquisition, supervision, independent
implementation of the real-time simulation, data analysis, writing

## References

  * (1) Brambilla, N., et al.: Heavy Quarkonium: Progress, Puzzles, and Opportunities. Eur. Phys. J. C71, 1534 (2011). doi:10.1140/epjc/s10052-010-1534-9. 1010.5827
  * (2) Rothkopf, A.: Heavy Quarkonium in Extreme Conditions. Phys. Rept. 858, 1–117 (2020). doi:10.1016/j.physrep.2020.02.006. 1912.02253
  * (3) Strickland, Michael: Pseudothermalization of the quark-gluon plasma. J. Phys. Conf. Ser. 1602(1), 012018 (2020). doi:10.1088/1742-6596/1602/1/012018. 2006.09284
  * (4) Berges, J., Heller, M.P., Mazeliauskas, A., Venugopalan, R.: Thermalization in QCD: theoretical approaches, phenomenological applications, and interdisciplinary connections (2020). 2005.12299
  * (5) Shen, C.: Studying QGP with flow: A theory overview. In: 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions (2020). 2001.11858
  * (6) Akamatsu, Y.: Approach to thermalization and hydrodynamics. In: 28th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions (2020). 2001.01429
  * (7) Kurkela, A., Mazeliauskas, A., Paquet, J.-F., Schlichting, S., Teaney, D.: Matching the non-equilibrium initial stage of heavy ion collisions to hydrodynamics with QCD kinetic theory. PoS Confinement2018, 152 (2019). doi:10.22323/1.336.0152
  * (8) Larsen, R., Meinel, S., Mukherjee, S., Petreczky, P.: Excited bottomonia in quark-gluon plasma from lattice QCD. Phys. Lett. B 800, 135119 (2020). doi:10.1016/j.physletb.2019.135119. 1910.07374
  * (9) Offler, S., Aarts, G., Allton, C., Glesaaen, J., Jäger, B., Kim, S., Lombardo, M.P., Ryan, S.M., Skullerud, J.-I.: News from bottomonium spectral functions in thermal QCD. PoS LATTICE2019, 076 (2019). doi:10.22323/1.363.0076. 1912.12900
  * (10) Larsen, R., Meinel, S., Mukherjee, S., Petreczky, P.: Thermal broadening of bottomonia: Lattice nonrelativistic QCD with extended operators. Phys. Rev. D 100(7), 074506 (2019). doi:10.1103/PhysRevD.100.074506. 1908.08437
  * (11) Kim, S., Petreczky, P., Rothkopf, A.: Quarkonium in-medium properties from realistic lattice NRQCD. JHEP 11, 088 (2018). doi:10.1007/JHEP11(2018)088. 1808.08781
  * (12) Burnier, Y., Ding, H.-T., Kaczmarek, O., Kruse, A.-L., Laine, M., Ohno, H., Sandmeyer, H.: Thermal quarkonium physics in the pseudoscalar channel. JHEP 11, 206 (2017). doi:10.1007/JHEP11(2017)206. 1709.07612
  * (13) Burnier, Y., Kaczmarek, O., Rothkopf, A.: Quarkonium at finite temperature: Towards realistic phenomenology from first principles. JHEP 12, 101 (2015). doi:10.1007/JHEP12(2015)101. 1509.07366
  * (14) Burnier, Y., Kaczmarek, O., Rothkopf, A.: In-medium P-wave quarkonium from the complex lattice QCD potential. JHEP 10, 032 (2016). doi:10.1007/JHEP10(2016)032. 1606.06211
  * (15) Song, T., Gubler, P., Hong, J., Lee, S.H., Morita, K.: $J/\psi$ near $T_{c}$ (2020). 2009.08741
  * (16) Braga, N.R.F., Ferreira, L.F., Vega, A.: Holographic model for charmonium dissociation. Phys. Lett. B774, 476–481 (2017). doi:10.1016/j.physletb.2017.10.013. 1709.05326
  * (17) Grigoryan, H.R., Hohler, P.M., Stephanov, M.A.: Towards the Gravity Dual of Quarkonium in the Strongly Coupled QCD Plasma. Phys. Rev. D82, 026005 (2010). doi:10.1103/PhysRevD.82.026005. 1003.1138
  * (18) Fujita, M., Kikuchi, T., Fukushima, K., Misumi, T., Murata, M.: Melting Spectral Functions of the Scalar and Vector Mesons in a Holographic QCD Model. Phys. Rev. D81, 065024 (2010). doi:10.1103/PhysRevD.81.065024. 0911.2298
  * (19) Borghini, N., Gombeaud, C.: Heavy quarkonia in a medium as a quantum dissipative system: Master equation approach. Eur. Phys. J. C72, 2000 (2012). doi:10.1140/epjc/s10052-012-2000-7. 1109.4271
  * (20) Akamatsu, Y., Rothkopf, A.: Stochastic potential and quantum decoherence of heavy quarkonium in the quark-gluon plasma. Phys. Rev. D85, 105011 (2012). doi:10.1103/PhysRevD.85.105011. 1110.1203
  * (21) Akamatsu, Y.: Real-time quantum dynamics of heavy quark systems at high temperature. Phys. Rev. D87(4), 045016 (2013). doi:10.1103/PhysRevD.87.045016. 1209.5068
  * (22) Brambilla, N., Escobedo, M.A., Strickland, M., Vairo, A., Vander Griend, P., Weber, J.H.: Bottomonium suppression in an open quantum system using the quantum trajectories method (2020). 2012.01240
  * (23) De Jong, W.A., Metcalf, M., Mulligan, J., Płoskoń, M., Ringer, F., Yao, X.: Quantum simulation of open quantum systems in heavy-ion collisions (2020). 2010.03571
  * (24) Yao, X., Mehen, T.: Quarkonium Semiclassical Transport in Quark-Gluon Plasma: Factorization and Quantum Correction (2020). 2009.02408
  * (25) Yao, X., Ke, W., Xu, Y., Bass, S.A., Müller, B.: Coupled Boltzmann Transport Equations of Heavy Quarks and Quarkonia in Quark-Gluon Plasma (2020). 2004.06746
  * (26) Ålund, O., Akamatsu, Y., Laurén, F., Miura, T., Nordström, J., Rothkopf, A.: Trace preserving quantum dynamics using a novel reparametrization-neutral summation-by-parts difference operator (2020). 2004.04406
  * (27) Miura, T., Akamatsu, Y., Asakawa, M., Rothkopf, A.: Quantum Brownian motion of a heavy quark pair in the quark-gluon plasma. Phys. Rev. D 101(3), 034011 (2020). doi:10.1103/PhysRevD.101.034011. 1908.06293
  * (28) Brambilla, N., Escobedo, M.A., Vairo, A., Vander Griend, P.: Transport coefficients from in medium quarkonium dynamics. Phys. Rev. D 100(5), 054025 (2019). doi:10.1103/PhysRevD.100.054025. 1903.08063
  * (29) Yao, X., Müller, B.: Quarkonium inside the quark-gluon plasma: Diffusion, dissociation, recombination, and energy loss. Phys. Rev. D100(1), 014008 (2019). doi:10.1103/PhysRevD.100.014008. 1811.09644
  * (30) Blaizot, J.-P., Escobedo, M.A.: Approach to equilibrium of a quarkonium in a quark-gluon plasma. Phys. Rev. D98(7), 074007 (2018). doi:10.1103/PhysRevD.98.074007. 1803.07996
  * (31) Brambilla, N., Escobedo, M.A., Soto, J., Vairo, A.: Heavy quarkonium suppression in a fireball. Phys. Rev. D 97(7), 074009 (2018). doi:10.1103/PhysRevD.97.074009. 1711.04515
  * (32) Kajimoto, S., Akamatsu, Y., Asakawa, M., Rothkopf, A.: Dynamical dissociation of quarkonia by wave function decoherence. Phys. Rev. D 97(1), 014003 (2018). doi:10.1103/PhysRevD.97.014003. 1705.03365
  * (33) Yao, X., Müller, B.: Approach to equilibrium of quarkonium in quark-gluon plasma. Phys. Rev. C97(1), 014908 (2018). doi:10.1103/PhysRevC.97.049903, 10.1103/PhysRevC.97.014908. [Erratum: Phys. Rev.C97,no.4,049903(2018)]. 1709.03529
  * (34) Akamatsu, Y.: Quarkonium in Quark-Gluon Plasma: Open Quantum System Approaches Re-examined (2020). 2009.10559
  * (35) Brambilla, N., Pineda, A., Soto, J., Vairo, A.: Effective field theories for heavy quarkonium. Rev. Mod. Phys. 77, 1423 (2005). doi:10.1103/RevModPhys.77.1423. hep-ph/0410047
  * (36) Laine, M., Philipsen, O., Romatschke, P., Tassler, M.: Real-time static potential in hot QCD. JHEP 03, 054 (2007). doi:10.1088/1126-6708/2007/03/054. hep-ph/0611300
  * (37) Beraudo, A., Blaizot, J.-P., Ratti, C.: Real and imaginary-time Q anti-Q correlators in a thermal medium. Nucl. Phys. A806, 312–338 (2008). doi:10.1016/j.nuclphysa.2008.03.001. 0712.4394
  * (38) Brambilla, N., Ghiglieri, J., Vairo, A., Petreczky, P.: Static quark-antiquark pairs at finite temperature. Phys. Rev. D78, 014017 (2008). doi:10.1103/PhysRevD.78.014017. 0804.0993
  * (39) Rothkopf, A., Hatsuda, T., Sasaki, S.: Complex Heavy-Quark Potential at Finite Temperature from Lattice QCD. Phys. Rev. Lett. 108, 162001 (2012). doi:10.1103/PhysRevLett.108.162001. 1108.1579
  * (40) Burnier, Y., Kaczmarek, O., Rothkopf, A.: Static quark-antiquark potential in the quark-gluon plasma from lattice QCD. Phys. Rev. Lett. 114(8), 082001 (2015). doi:10.1103/PhysRevLett.114.082001. 1410.2546
  * (41) Burnier, Y., Rothkopf, A.: Complex heavy-quark potential and Debye mass in a gluonic medium from lattice QCD. Phys. Rev. D 95(5), 054511 (2017). doi:10.1103/PhysRevD.95.054511. 1607.04049
  * (42) Petreczky, P., Weber, J.: Lattice Calculations of Heavy Quark Potential at Finite Temperature. Nucl. Phys. A 967, 592–595 (2017). doi:10.1016/j.nuclphysa.2017.04.011. 1704.08573
  * (43) Petreczky, P., Rothkopf, A., Weber, J.: Realistic in-medium heavy-quark potential from high statistics lattice QCD simulations. Nucl. Phys. A 982, 735–738 (2019). doi:10.1016/j.nuclphysa.2018.10.012. 1810.02230
  * (44) Lafferty, D., Rothkopf, A.: Improved Gauss law model and in-medium heavy quarkonium at finite density and velocity. Phys. Rev. D 101(5), 056010 (2020). doi:10.1103/PhysRevD.101.056010. 1906.00035
  * (45) Burnier, Y., Rothkopf, A.: Disentangling the timescales behind the non-perturbative heavy quark potential. Phys. Rev. D 86, 051503 (2012). doi:10.1103/PhysRevD.86.051503. 1208.1899
  * (46) Burnier, Y., Rothkopf, A.: A hard thermal loop benchmark for the extraction of the nonperturbative $Q\bar{Q}$ potential. Phys. Rev. D 87, 114019 (2013). doi:10.1103/PhysRevD.87.114019. 1304.4154
  * (47) Rothkopf, A.: Bayesian techniques and applications to QCD. PoS Confinement2018, 026 (2018). doi:10.22323/1.336.0026. 1903.02293
  * (48) D’Onofrio, M., Rummukainen, K., Tranberg, A.: Sphaleron Rate in the Minimal Standard Model. Phys. Rev. Lett. 113(14), 141602 (2014). doi:10.1103/PhysRevLett.113.141602. 1404.3565
  * (49) Micha, R., Tkachev, I.I.: Turbulent thermalization. Phys. Rev. D 70, 043538 (2004). doi:10.1103/PhysRevD.70.043538. hep-ph/0403101
  * (50) Berges, J., Rothkopf, A., Schmidt, J.: Non-thermal fixed points: Effective weak-coupling for strongly correlated systems far from equilibrium. Phys. Rev. Lett. 101, 041603 (2008). doi:10.1103/PhysRevLett.101.041603. 0803.0131
  * (51) Berges, J., Wallisch, B.: Nonthermal Fixed Points in Quantum Field Theory Beyond the Weak-Coupling Limit. Phys. Rev. D 95(3), 036016 (2017). doi:10.1103/PhysRevD.95.036016. 1607.02160
  * (52) Chatrchyan, A., Geier, K.T., Oberthaler, M.K., Berges, J., Hauke, P.: Analog reheating of the early universe in the laboratory (2020). 2008.02290
  * (53) Laine, M., Philipsen, O., Tassler, M.: Thermal imaginary part of a real-time static potential from classical lattice gauge theory simulations. JHEP 09, 066 (2007). doi:10.1088/1126-6708/2007/09/066. 0707.2458
  * (54) Lehmann, A.: Minkowskian lattice simulation for non-relativistic quarks in classical fields. PhD thesis, Heidelberg University, Heidelberg (June 2020). doi:10.11588/heidok.00028377
  * (55) Arnold, P.: Hotbviolation, the lattice, and hard thermal loops. Physical Review D 55(12), 7781–7796 (1997). doi:10.1103/physrevd.55.7781
  * (56) WATSON, G.N.: Three triple integrals. The Quarterly Journal of Mathematics os-10(1), 266–276 (1939). doi:10.1093/qmath/os-10.1.266
  * (57) Glasser, M.L., Boersma, J.: Exact values for the cubic lattice green functions. Journal of Physics A: Mathematical and General 33(28), 5017–5023 (2000). doi:10.1088/0305-4470/33/28/306
  * (58) Bodeker, D.: On the effective dynamics of soft nonAbelian gauge fields at finite temperature. Phys. Lett. B 426, 351–360 (1998). doi:10.1016/S0370-2693(98)00279-2. hep-ph/9801430
  * (59) Bodeker, D., Moore, G.D., Rummukainen, K.: Chern-Simons number diffusion and hard thermal loops on the lattice. Phys. Rev. D 61, 056003 (2000). doi:10.1103/PhysRevD.61.056003. hep-ph/9907545
  * (60) Hückel, E., Debye, P.: The theory of electrolytes: I. lowering of freezing point and related phenomena. Phys. Z 24, 185–206 (1923)
  * (61) Kasper, V., Hebenstreit, F., Berges, J.: Fermion production from real-time lattice gauge theory in the classical-statistical regime. Physical Review D 90(2), 025016 (2014). doi:10.1103/PhysRevD.90.025016. arXiv: 1403.4849. Accessed 2020-12-06
  * (62) Wilson, K.G.: Confinement of quarks. Physical Review D 10(8), 2445–2459 (1974). doi:10.1103/physrevd.10.2445
  * (63) Klassen, T.R.: The anisotropic wilson gauge action. Nuclear Physics B 533(1-3), 557–575 (1998). doi:10.1016/s0550-3213(98)00510-0
  * (64) Ambjorn, J., Askgaard, T., Porter, H., Shaposhnikov, M.E.: Sphaleron transitions and baryon asymmetry: A Numerical real time analysis. Nucl. Phys. B 353, 346–378 (1991). doi:10.1016/0550-3213(91)90341-T
  * (65) Ambjorn, J., Krasnitz, A.: Improved determination of the classical sphaleron transition rate. Nucl. Phys. B 506, 387–403 (1997). doi:10.1016/S0550-3213(97)00524-5. hep-ph/9705380
  * (66) Grigoriev, D.Y., Rubakov, V.A.: Soliton Pair Creation at Finite Temperatures. Numerical Study in (1+1)-dimensions. Nucl. Phys. B 299, 67–78 (1988). doi:10.1016/0550-3213(88)90466-X
  * (67) Akamatsu, Y., Rothkopf, A., Yamamoto, N.: Non-Abelian chiral instabilities at high temperature on the lattice. JHEP 03, 210 (2016). doi:10.1007/JHEP03(2016)210. 1512.02374
  * (68) Blackman, R.B., Tukey, J.W.: The measurement of power spectra from the point of view of communications engineering — part i. The Bell System Technical Journal 37(1), 185–282 (1958). doi:10.1002/j.1538-7305.1958.tb03874.x
  * (69) Rothkopf, A.: ”Classical Statistical Simulation of SU(3) Yang-Mills Theory in the Presence of Static Sources”. doi:10.5281/zenodo.4332406. https://doi.org/10.5281/zenodo.4332406
  * (70) Yanagihara, R., Iritani, T., Kitazawa, M., Asakawa, M., Hatsuda, T.: Distribution of Stress Tensor around Static Quark–Anti-Quark from Yang-Mills Gradient Flow. Physics Letters B 789, 210–214 (2019). doi:10.1016/j.physletb.2018.09.067. arXiv: 1803.05656. Accessed 2020-03-21
  * (71) Lehmann, A., Rothkopf, A.: Real-Time-Evolution of Heavy-Quarkonium Bound States. PoS LATTICE2019, 074 (2019). doi:10.22323/1.363.0074. 2003.02509