11institutetext: Faculty of Science and Technology, University of Stavanger,
NO-4021 Stavanger, Norway

# Inverse problems, real-time dynamics and lattice simulations

Alexander Rothkopf 1133<EMAIL_ADDRESS>
###### Abstract

The determination of real-time dynamics of strongly coupled quantum fields is
a central goal of modern nuclear and particle physics, which requires insight
into quantum field theory beyond the weak-coupling approximation. While
lattice QCD has provided vital insights into the non-perturbative static
properties of quarks and gluons it hides their real-time dynamics behind an
ill-posed inverse problem. In this proceeding I will discuss developments in
tackling the inverse problem on the lattice and touch upon progress in the
direct simualtion of real-time dynamics.

##  1 Introduction

The task of determining the real-time dynamics of strongly correlated quantum
fields arises in many of the topical tracks of this conference, be it in the
study of the transport properties of the hot quark-gluon plasma created in
relativistic heavy-ion collisions in Track D or the relaxation of the dense
interior of collisions of two neutron stars in Track F. With the electron-ion
collider project underway the phase-space tomography of nucleons and
eventually nuclei constitutes a central topic in Track A and Track B.
Similarly the precision determination of particle properties, such as the muon
$g-2$ requires high precision real-time input from the theory of the strong
interactions, quantum chromodynamics (QCD), as discussed in contributions to
Track E.

Lattice QCD over the past decade has proven to be a robust tool to elucidate
strong interaction physics (see e.g. Ratti:2018ksb  or Davoudi:2022bnl ). One
sets out to evaluate expectation values of observables in Feynman’s path
integral from first-principles on a hypercubic grid with lattice spacing $a$

$\displaystyle\langle{\cal O}(x_{1}){\cal O}(x_{2})\rangle=\int{\cal
D}[U,\bar{\psi},\psi]\,{\cal O}[x_{1};U,\bar{\psi},\psi]{\cal
O}[x_{2};U,\bar{\psi},\psi]\,{\rm exp}\Big{[}iS_{\rm
QCD}[U,\bar{\psi},\psi]\Big{]}.$ (1)

In order to preserve the gauge invariance of the theory after discretization,
the $su(3)$ valued gauge fields are formulated in terms of so called link
variables $U_{\mu}(x)={\rm exp}[-igA^{\mu}(x+a\hat{\mu}/2)]$ located in
between the nodes of the grid. Quark fields $\psi,\bar{\psi}$ on the other
hand are distributed on the nodes.

In its Minkowski time form the Feynman weight ${\rm exp}[-iS]$ amounts to a
pure phase, rendering the evaluation of this path integral by means of naive
statistical sampling exponentially hard. The success underlying lattice QCD
simulations is the realization that we may analytically continue the real-time
path integral onto the negative imaginary time axis, i.e. into the Euclidean
domain. There the Feynman weight ${\rm exp}[-S_{E}]$ is real-valued and
bounded from below, which offers a statistical interpretation of the path
integral and makes its evaluation amenable to established sampling approaches,
such as Markov-Chain Monte-Carlo. I.e one generates a collection of field
realizations distributed according to the Feynman weight, computes the value
of operators on each of these, and by taking the mean systematically
approximates the quantum-statistical expectation value of interest. Due to the
subtle but vital relation between the physical extent of the Euclidean domain
and the inverse temperature of the system, introducing periodic boundary
condition in imaginary time is not a discretization artifact but actually
encodes a vital physical property of the system.

Our original intention is to evaluate Minkowski-time correlation functions,
which in Euclidean time on the lattice are not directly accessible. One
central player in the quest of extracting such real-time information is the
spectral function. The spectral function $\rho(\omega)$ of a two-point
function refers to the unique quantity, which, through an integral transform
with a known kernel function $K$, encodes all of the information of the
correlator in a representation independent fashion, i.e. irrespective of
whether the two-point function is expressed as retarded or advanced real-time
quantity or on the Euclidean domain

$\displaystyle{\cal O}^{E}(\tau)=\int
d\omega\,K(\omega,\tau)\rho(\omega)\quad\overset{\rm
analytic\,continuation}{\leftrightarrow}\quad\int d\omega\,K(\omega,it\pm
i\varepsilon)\rho(\omega)={\cal O}^{R/A}(t).$ (2)

Importantly the spectral function itself often encodes in an intuitive fashion
the physics of interest. Take as an example the heavy quarkonium spectral
functions studied in Track C and Track D. The low frequency behavior of
$\rho(\omega)$ encodes transport properties, such as heavy quark diffusion via
linear response theory (Kubo formula) and at intermediate frequencies the in-
medium hadron bound-state properties are encoded in peak structures (for a
more detailed discussion see e.g. Rothkopf:2019ipj ).

##  2 Spectral function reconstruction from Euclidean lattices

We thus face the task of extracting spectral functions from Euclidean
simulations. It turns out that lattice QCD acts on spectral functions very
similarly to an imperfect detector acts on particle tracks. The integral
transform of eq. 2 smears out the intricate structures present in the spectral
function (see fig. 1) into an imaginary time correlator $D$, based on which we
must solve an inverse problem to regain access to the values of $\rho$. The
challenge lies in the fact that a Monte-Carlo simulation produces only
discrete and noisy realizations of the correlation functions and that the
kernel in the convolution can lead to exponential loss of information when
going from $\rho$ to $D$.

Figure 1: Sketch of the inverse problem using quarkonium as example. We must
extract spectral functions $\rho$ from Euclidean lattice simulations, which
only produce a noisy and discrete imaginary time correlator $D$ that
corresponds to a convolution of the spectral function with a known kernel
function K.

In some circumstances the inverse problem to reconstruct $\rho$ is benign,
when the spectral structures are well separated and their form is known
apriori. In that case, e.g. when extracting the grounds state masses of
hadrons at $T=0$, a parametrized model of the spectral structure can be
constructed which is directly fitted to the correlator $D$. In other
circumstances, e.g. when extracting in-medium hadron properties, intricate
features appear in the spectral function that are not well separated and no
simple closed form parametrization of $\rho$ is available. In that case the
inverse problem is clearly ill-posed in the sense of Hadamard: there exist
many different realizations of $\rho$ which reproduce the noisy correlator $D$
within its statistical uncertainties (non-uniqueness) and due to the
exponential behavior of the kernel $K$ the inverse problem without additional
regularization is highly sensitive to noise in the input data (instability).

The exponential information loss induced by the integral kernel $K$ is a
common issue when reconstructing spectral functions in different contexts. It
is obvious from the functional form of $K$ when considering hadron physics at
$T=0$ and $T>0$

$\displaystyle K^{T=0}={\rm exp}[-\omega\tau],\quad K^{T>0}={\rm
cosh}[\omega(\tau-\beta/2)]/{\rm sinh}[\omega\beta/2],\quad K^{\rm PDF}_{\rm
Ioffe}={\rm cos}[x\nu]$ (3)

However one faces an additional challenge, due to limited reach of the
simulations. At finite temperature, we are apriori limited to accessing
Euclidean times only up to $\tau<\beta$, which leads to the inherent
coarseness of Matsubara frequencies. It is this coarseness which may prevent
the successful reconstruction of essential spectral features related to
thermal physics. In the study of PDFs in the pseudo-PDF formalism
HadStruc:2021qdf  one is in practice limited to resolving only a subset of
Ioffe-time values, which leads to exponentially small eigenvalues in the
kernel matrix $K$ associated with the otherwise innocent looking kernel
${\cos}[\nu x]$. Keeping these challenges in mind, we must ask how to extract
most accurately the information residing inside the simulated correlators.

Three strategies can be found in the literature: Strategy I does not accept
the premise that lattice QCD can provide us with "the spectral function".
Instead it settles for extracting a smeared approximation of the spectral
function from a linear combination of the input data. This strategy is
implemented using the modified Backus-Gilbert method backus1968resolving ;
hansen2019extraction . Strategy II accepts the premise and attempts to model
the spectral structures using prior domain information input, derived e.g.
from low-energy effective field theories or perturbation theory. This leads to
a sophisticated ansatz for the spectral function fitted to the correlator (see
e.g. mages2015shear ; burnier2017thermal ; Astrakhantsev:2018oue ). Strategy
III also accepts the premise but refrains from a direct parametrization. It
includes prior information on the level of a regulator functional introduced
in the context of Bayes theorem. The Bayesian inference approach encompasses
various methods from sparse modelling, Tikhonov regularization
tikhonov1943stability  and even some neural network approaches.

In the following I will focus on Bayesian inference jarrell1996bayesian , a
framework which is also applied in various contexts in the contributions to
Track H. The discrete spectral relation we wish to invert reads

$\displaystyle
D_{i}=\sum_{k=1}^{N_{\omega}}\,\Delta\omega_{k}K(\omega_{k},\tau_{i})\,\rho(\omega_{i})\quad
i\in[1,\ldots,N_{\tau}],\;N_{\omega}\gg N_{\tau}.$ (4)

I.e. we wish to reconstruct the $N_{\omega}$ parameters
$\rho_{k}=\rho(\omega_{k})$ from $N_{\tau}$ noisy datapoints. To give meaning
to this task, one deploys Bayes theorem, which describes the conditional
probability for a test function $\rho$ to be the correct spectral function
underlying supplied simulation data $D$ and prior information $I$. This so
called posterior probability is rewritten as a product

$\displaystyle\underbracket{P[\rho|D,I]}_{\rm
posterior}=\underbracket{P[D,\rho,I]}_{\rm
likelihood}\underbracket{P[\rho|I]}_{\rm prior}/\underbracket{P[D|I]}_{\rm
evidence}\quad\overset{\rm MAP}{\rightarrow}\quad\left.\delta
P[\rho|D,I]/\delta\rho\right|_{\rho=\rho^{\rm MAP}}=0$ (5)

of the likelihood probability $P[D,\rho,I]\propto{\rm exp}[-L]$, which encodes
information about the data generation process and the prior probability
$P[\rho|I]\propto{\rm exp}[S]$, which can be used to encode domain knowledge
beyond what is provided by the simulation data. Normalization is provided by
the so-called evidence. In order to save computational cost, one often does
not sample to full posterior but instead computes the most likely spectral
function $\rho^{\rm MAP}$, the maximum aposteriori, which denotes the extremum
of the posterior. In the absence of prior information one ends up with many
degenerate maximum likelihood solutions. The functional $S$ inside the prior
is therefore called regulator, as it picks from among the many degenerate
local extrema of the likelihood the one that aligns closest with given domain
information (e.g. positivity of a hadronic spectral function).

The prior probability is conventionally characterized by its two first
cumulants, the so-called default model $m_{l}$ and (hyperparameter) weight
$\alpha_{l}$. Whereas $m_{l}$ encodes the most probable apriori value of the
of the parameter $\rho_{l}$ its apriori variance is described by
$1/\alpha_{l}$. All Bayesian approaches by construction must agree in the
Bayesian continuum limit of simultaneously increasing the number of data
points $N_{\tau}\to\infty$ and reducing their uncertainty $\Delta D/D\to 0$.
The choice of $S$ influences which types of artifacts one encounters on the
way, such as e.g. ringing for weak regulators or overdamping of spectral
structures for strong regulators.

Its versatility allows the Bayesian approach to encompass several
regularization schemes. The well-known Tikhonov approach is found to
correspond to a Gaussian prior probability with a vanishing default model and
a weight $\alpha$ set to a value that recovers $L=N_{\tau}/2$ for $\rho^{\rm
MAP}$. Note that this approach does not exploit the positivity of hadronic
spectral functions as its regulators admits also negative values for the
$\rho_{l}$’s. Two other Bayesian approaches, which are commonly deployed in
the study of positive definite spectral functions are the Maximum Entropy
Method (MEM) jarrell1996bayesian ; asakawa2001maximum  and the Bayesian
Reconstruction (BR) burnier2013bayesian  method, which feature the following
regulators

$\displaystyle S_{MEM}=\alpha\int d\omega\Big{(}\rho-m-\rho{\rm
log}\Big{[}\frac{\rho}{m}\Big{]}\Big{)},\quad S_{BR}=\alpha\int
d\omega\Big{(}1-\frac{\rho}{m}+{\rm log}\Big{[}\frac{\rho}{m}\Big{]}\Big{)}.$
(6)

The former was originally constructed for the reconstruction of images in
astronomy and is based on four axioms. For a detailed discussion of the axioms
see skilling1988axioms . The hyperparameter $\alpha$ is treated in a self
consistent fashion in the literature by finding the extremum of the evidence,
a procedure that relies on a saddle-point approximation of the posterior,
which is otherwise analytically intractable.

The BR method is also based on four axioms burnier2013bayesian , which are
selected with the one-dimensional reconstruction problem in lattice QCD in
mind. One of the axioms is scale invariance, which in essence requires that
the end result of the reconstruction may not depend on the units of the
spectral function and is takes care of the fact that $\rho$ is not necessarily
a probability distribution itself. The second axiom to note is smoothness,
which requires that the end result should be a smooth (at least twice
differentiable) function. The functional form of $S_{BR}$ allows for a semi-
analytic apriori marginalization of $\alpha$, which offers a truly Bayesian
treatment of this often poorly constrained hyperparameter.

Bayesian spectral function reconstruction has been deployed to the study of
$T>0$ physics and $T=0$ physics from the lattice. The study of in-medium
hadron spectral functions, in particular those of heavy-quarkonium bound
states Aarts:2014cda ; Kim:2018yhk  is a central application, with new results
presented by the FASTSUM collaboration in Track D at this conference (for a
recent study of excited states quarkonium see Larsen:2019zqv ). Progress has
been made by using Bayesian spectral reconstructions in the study of the
complex in-medium potential between static quarks from lattice QCD
Burnier:2015tda ; Burnier:2016mxc . Recent results on the potential from
lattices with dynamical quarks Bala:2021fkm  will be presented by Gaurang
Parkar in Track C. Over the past years first steps have been taken to extract
parton distribution functions from the lattice using among others Bayesian
inference Liang:2019frk ; Karpie:2019eiq , an overview over the state-of-the-
art in PDFs from the lattice is presented in Track B.

The past years have also seen heightened interest in machine learning
approaches to tackle the spectral function reconstruction challenge. In the
literature one finds direct neural network approaches, where methods used in
image reconstruction are transferred to spectral functions. The situation here
is often actually simpler than in image reconstruction since the decoder step
of going from the spectral function to the data is analytically known. Two
recent approaches are supervised Kernel-Ridge regression Offler:2021fmg  or
support vector machine regression Kades:2019wtd ; fournier2020artificial ,
which differ by the basis functions used, the loss functional they deploy and
the regulators included to obtain unique solution.

A different approach to regularization has been explored in recent papers
Karpie:2019eiq ; Shi:2022yqw  where a neural network is used in the
parametrization of the spectral function itself $\rho(\omega)={\rm
NN}(\omega)$. This parametrization, which in essence amounts to a certain
choice of basis functions can be used in both supervised and unsupervised
settings. It is important to note that the regularization of the
reconstruction is still implemented through additional terms in the loss
functional that guides the training of these neural networks.

A very recent study Horak:2021syv  has turned to another well established
machine-learning technique, Gaussian processes, to reconstruct spectral
functions. It models both observed data and predicted data as arising from the
same underlying Gaussian distribution. It is the correlation matrix of this
high dimensional normal distribution that implements the regularization of the
approach. The hyperparameters of the correlation matrix are self-consistently
tuned based on training data.

If we take a step back and ask, where do we stand in regard to spectral
function reconstruction we can recognize that so far no clear methods winner
has emerged. The central limitation remains the scarce information content of
the correlator and one pressing open question is: how do we include more QCD
specific prior knowledge into the regularization of the inversion task? It is
important to note that significant increases in statistics by now enable the
use of complimentary methods, such as e.g. the Padé reconstruction
Tripolt:2018xeo .

As the community works on improving lattice simulations to more efficiently
approach the continuum limit (Szymanzik improvement) one encounters competing
interests: while improved actions may accelerate the approach to the
continuum, they can interfere with conventional spectral decomposition by
introducing e.g. negative spectral weights (see discussion in Bala:2021fkm ).
A better understanding of these artifacts may allow us to incorporate and
potentially correct for them in improved regulator terms.

Given the challenges faced in spectral function reconstruction from Euclidean
lattices, one may ask whether direct simulation strategies for real-time
dynamics can provide complementary insight.

##  3 Direct real-time simulations

The notorious sign problem gattringer_approaches_2016  prevents a straight
forward evaluation of the Minkowski-time path integral using established
Monte-Carlo approaches. On the one hand we may find scenarios, where the
quantum fluctuations of fields are small compared to the statistical
fluctuations. This approximation underlies the classical statistical lattice
gauge theory (CLGT) framework, which has been successfully used to shed light
on the non-equilibrium physics of the early stages of heavy-ion collisions
(see e.g. Berges:2013fga ), as well as on topology changing processes in
thermal equilibrium (see e.g. Moore:2010jd ). In a classical statistical
simulation, the gauge fields are evolved deterministically according to the
classical equations of motion, initialized by values drawn from a thermal
distribution. Carrying out the deterministic evolution for multiple different
stochastic initial values allows us to form ensemble averages of observables
at finite real-time.

On the other hand one may attempt to design alternative techniques for the
direct evaluation of the real-time path integral. The two most promising
strategies in our community are the Lefschetz thimbles cristoforetti_new_2012
and complex Langevin namiki_stochastic_1992 , both of which rely on a
complexification of the field degrees of freedom, moving the path integration
into the complex plane. The challenge with Lefschetz thimbles lies in
identifying the corresponding submanifolds, by solving a flow equation, which
can incur large numerical cost in practice. For complex Langevin, the main
drawback lies in the fact that while the approach may converge, it does not
necessarily converge to the correct solution.

Taking up a line of inquiry started in Laine:2007qy , we recently used
classical statistical simulations to explore the binding of static quark-
antiquark pairs in real-time Lehmann:2020fjt . For the first time it was
possible to compute the real part of the interaction potential by
incorporating the physics of the sources into the Gauss law of the simulation.
The static potential is related to the time evolution of the Wilson loop,
which describes the propagation of a pair of a static quark and antiquark. It
has been shown on general grounds Burnier:2012az  that if the Wilson loop
evolves according to a Schrödinger equation at late times, its spectral
function, which in Minkowski time is obtained from a simple Fourier transform
of the real-time Wilson loop, will exhibit a dominant skewed Lorentzian peak

$\displaystyle i\partial_{t}\langle W_{\square}(r,t)\rangle\sim\Big{(}{\rm
Re}$ $\displaystyle[V](r)-i|{\rm Im}[V](r)|\Big{)}\langle
W_{\square}(r,t)\rangle,\qquad\overset{\rm Fourier}{\leftrightarrow}$
$\displaystyle\langle W_{\square}(r,\omega)\rangle\sim c_{0}\frac{{\rm
Im}[V]^{2}(r)}{(\omega-{\rm Re}[V](r))^{2}+{\rm Im}[V]^{2}(r)}+{\rm
background}.$ (7)

I.e. the position and width of the spectral peak encodes the real- and
imaginary part of the complex potential. Note that the evaluation of the
Wilson loop in a fully quantum path integral or as an observable in a
classical statistical simulation differs substantially. In the former setting
the evaluation of the Wilson loop actually corresponds to a reweighting of a
system without static quarks to a system where the sources are present, i.e.
$\langle W(r,\tau)\rangle_{\rm path\,integral}=Z_{\rm
medium\,+\,Q\bar{Q}}(r,\tau)/Z_{\rm medium}$. in classical statistical
simulations on the other hand the Wilson loop is simply evaluated as
observable and does not affect the real-time dynamics of the gauge fields.
Sources instead need to be treated as explicit part of the equations of
motion, where their charge density $\sigma$ enters as

$\displaystyle Z_{\rm CLGT}$ $\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{)},$
(8) $\displaystyle j^{0}$ $\displaystyle=M\Big{(}\delta^{(3)}({\bf x}-{\bf
x}_{0})-\delta^{(3)}({\bf x}-{\bf
x}_{1})\Big{)}={\color[rgb]{1,0,0}\sigma({\bf x};{\bf x}_{0},{\bf x}_{1})}.$
(9)

I.e. we must modify the local Gauss law on the lattice by hand, when
introducing static sources in the classical statistical simulation

$\displaystyle G({\bf x},t)\equiv\sum_{i}\Big{[}E_{i}({\bf x},t)-U_{i}({\bf
x},t)E_{i}({\bf x}+{\bf i},t)U_{-i}^{\dagger}({\bf
x},t)\Big{]}={\color[rgb]{1,0,0}\sigma({\bf x};{\bf x}_{0},{\bf x}_{1})}.$
(10)

This version of Gauss’ law amended by the static charge density enters the
classical statistical dynamics solely via the initial conditions, which must
be projected onto the constraint surface given by $G({\bf
x},t)={\color[rgb]{1,0,0}\sigma({\bf x};{\bf x}_{0},{\bf x}_{1})}$.

Solving the Yang-Mills equations of motion in the presence of the proper Gauss
law significantly modifies the behavior of the Wilson loop in CLGT. In
previous studies, where the charge density was absent in Gauss’ law, the
Wilson loop remained purely real and exhibited a simple exponential falloff at
late real-times. As shown in the left panel of fig. 2 in the presence of the
proper Gauss law, the real part of $\langle W_{\square}(r,t)\rangle$ shows an
intricate oscillatory behavior with a continuously damped amplitude over
time.The different colored curves correspond to different spatial separation
distances, the further apart the sources the faster the damping.

Figure 2: (left) Real part of the real-time Wilson loop in CLGT in the
presence of the proper Gauss law. Note the both oscillatory and damped
component of the evolution. Its imaginary part shows a similar oscillatory
pattern. (right) The spectral function of the real-time Wilson loop, which
exhibits a well-defined low-lying peak of skewed-Lorentzian shape. Figures
from Lehmann:2020fjt .

Inspection of the imaginary part of $\langle W_{\square}(r,t)\rangle$ (see
Lehmann:2020fjt ) further reveals that it differs from zero and it too
exhibits damped oscillations. When one takes the (windowed) Fourier transform
of the Wilson loop at different spatial separation distances, the spectral
functions in the right panel of fig. 2 ensue. We see that the spectral
functions exhibit a well defined lowest-lying peak structure, whose peak
position and width change with the spatial separation of the static sources.
Interpreted in the context of the complex interquark potential this result
tells us that a finite real-part of the potential exists (finite position of
the peak), as well as a finite imaginary part is present (finite width of the
peak).

Figure 3: (left) Real part of the static interquark potential in thermal
classical lattice gauge theory Lehmann:2020fjt , which is compatible with the
expected Debye screening behavior in a gluonic plasma at high temperatures.
(right) A high precision analysis of the imaginary part of the interquark
potential in CLGT on large lattices Boguslavski:2021zga  and a comparison to a
fit motivated by analytic predictions from classical hard thermal loops.

Plotted against separation distance, we show in the left plot of fig. 3 the
extracted values of the real-part of the static potential in CLGT. This is the
first time that this quantity has been determined in classical statistical
lattice simulations. We find that its values are compatible with an
asymptotically flat form ${\rm Re}[V](r)\sim-{\rm exp}[-m_{D}r]/r+c$, which is
reminiscent of the expected Debye screening behavior for a gluon gas at high
temperature. On the left we show a recent high precision determination of the
imaginary part of the static potential from a study by a collaboration between
colleagues at Kent State University and TU Vienna Boguslavski:2021zga . The
large lattices deployed allow the authors to explore small enough distances to
match the simulations to continuum computations of the potential from
classical hard-thermal loops.

The relevance of a correct treatment of Gauss’ law on the lattice hints at
that care must be taken when designing improved actions applicable to systems
where translation invariance is broken e.g. by the presence of static charges.
In Rothkopf:2021jye  I discuss some of the insight the computational classical
electrodynamics community has gathered over the past decades. There it is
recognized that two concepts are central to establish accurate finite
difference based discretizations schemes. One is the use of so-called
summation by parts finite difference operators (see also Rothkopf:2022zfb ),
which mimic integration by parts in the discrete setting. In addition it is
known that a higher order discretization of Gauss’ law with a symmetric
central finite difference operators requires a so-called finite volume scheme,
which realizes not only the differential but also the integral form of Gauss’
law. Expressed in the language of finite differences the finite-volume scheme
amounts to the introduction of distributed sources, revealing an interesting
relation between the concept of smeared sources and the underlying
discretization of the field theory. In fig. 4 the failure of the naive central
finite difference scheme for the Abelian Gauss law is shown in the left panel,
where in the presence of a charge and anti-charge an artificial staggered
pattern of field line emerges (green arrows). Only after going over to a
finite volume scheme on the right, corresponding to a distributed charge, do
we obtain closed field lines (dark red) that accurately reproduce the correct
solution (light red arrows) and which satisfy the integral form of Gauss’ law
within the accuracy of the approximation.

Figure 4: Example of the difficulty encountered in discretizing Gauss’ law
with central finite differences. (left) Using a naive central finite
difference discretization an artificial staggered pattern of field lines
(green arrows) emerges which does not accurately reproduce the correct
solution (light red arrows). (right) A finite volume discretization based on
distributed sources achieves closed field lines (dark red arrows) and
conserves the integral form of Gauss’ law within the accuracy of the
discretization.

Last but not least let me discuss recent developments in the direct simulation
of real-time dynamics based on the complex Langevin framework
namiki_stochastic_1992 . One sets out to construct a stochastic process in an
artificial additional Langevin time dimension $\tau_{L}$ (not related to the
Euclidean time $\tau$), in which the real- and imaginary part of complexified
field degrees of freedom $\phi=\phi^{R}+i\phi^{I}$ are evolved in the presence
of Gaussian noise $\eta$

$\displaystyle\frac{d\phi^{R}}{d\tau_{L}}={\rm Re}\left[\left.i\frac{\delta
S[\phi]}{\delta\phi}\right|_{\phi=\phi^{R}+i\phi^{I}}\right]+\eta(\tau_{L}),\quad\frac{d\phi^{I}}{d\tau_{L}}={\rm
Im}\left[\left.i\frac{\delta
S[\phi]}{\delta\phi}\right|_{\phi=\phi^{R}+i\phi^{I}}\right].$ (11)

The goal is to produce in the late time limit a distribution
$P[\phi^{R},\phi^{I}]$ that, when used to evaluate expectation values,
reproduces the results of the Minkowski time path integral

$\displaystyle\int d\phi^{R}\int d\phi^{I}\,P[\phi^{R},\phi^{I}]\,{\cal
O}[\phi^{R}+i\phi^{I}]\leftrightarrow\int{\cal D}\phi\,{\cal O}[\phi]\,{\rm
exp}[iS].$ (12)

Traditionally there exist two major challenges: the occurrence of so-called
runaway trajectories, where the stochastic process spontaneously diverges and
the more fundamental limitation of convergence to incorrect results.

Recently my group has made progress with respect to runaways, by introducing
inherently stable implicit solvers to the simulation of complex Langevin
Alvestad:2021hsi . While previously runaways were combated by the introduction
of an adaptive step size into explicit solvers AartsJames2010 , the implicit
schemes offer additional benefits beyond inherent stability. It turns out that
they also provide an inherent regularization of the underlying path integral.
This has allowed us for the first time to simulate directly on the actual
real-time Schwinger-Keldysh contour without the need to introduce a tilt,
which so far has prevented access to phenomenologically relevant correlation
functions involving the degrees of freedom on the backward contour.

At this conference we will present our most recent work Alvestad:2022hsi  that
attacks the remaining challenge of complex Langevin, the possibility to
converge to the wrong solution. In Track G Daniel Alvestad will present our
novel strategy that exploits the freedom to introduce kernels into the complex
Langevin equation to restore correct convergence. Our approach acknowledges
that since the sign-problem is NP-hard Troyer:2004ge , no general solution
strategy is likely to succeed and therefore system specific prior knowledge
needs to be incorporated into the simulation. We do so by devising a learning
strategy that finds optimal kernels based on prior knowledge, such as the
symmetries of the system and correlation functions accessible in Euclidean
time. Our strategy allows us to extend threefold the real-time extend in which
complex Langevin simulations correctly reproduce the physics of the community
benchmark system, the strongly coupled anharmonic oscillator BergesSexty2007 .

The TU-Vienna group will also present work in progress Boguslavski:2022vjz  on
the simulation of real-time gauge theory on the lattice. They apply the gauge
cooling and dynamic stabilization strategy developed in the context of complex
Langevin for QCD at finite Baryo-chemical potential to simulations on the
tilted Schwinger-Keldysh contour. For a recent proposal on how to improve
real-time simulations based on the Lefschetz thimble approach see
Woodward:2022pet .

##  4 Conclusion

The determination of the real-time dynamics of strongly interacting matter
lies at the heart of many open questions discussed at this conference. Lattice
QCD simulations in Euclidean time have proven to be a robust tool for the
study of the static properties of the strong interactions but the ill-posed
inverse problem of extracting spectral functions severely limits our access to
real-time dynamics. The limited information content of the Euclidean
correlators presents a central challenge for conventional approaches such as
Bayesian inference or more recent machine learning strategies. Incorporation
of more specific QCD based prior information in the regularization of the
inverse problem is thus called for.

Direct real-time simulations offer an alternative route towards gaining
insight. Classical statistical lattice gauge theory neglects quantum
fluctuations and may provide some qualitative insight into e.g. the binding
properties of hadrons at very high temperatures. Complex Langevin simulations
are steadily improving and interest in their application to the real-time
setting has recently increased. Progress has been made to avoid runaways
apriori using implicit solvers and new strategies to overcome the convergence
to incorrect solution are actively explored.

##  5 Acknowledgments

A.R. is supported by the Research Council of Norway under the FRIPRO Young
Research Talent grant 286883.

## References

  * (1) C. Ratti, Rept. Prog. Phys. 81, 084301 (2018), 1804.07810
  * (2) Z. Davoudi et al., _Report of the Snowmass 2021 Topical Group on Lattice Gauge Theory_ , in _2022 Snowmass Summer Study_ (2022), 2209.10758
  * (3) A. Rothkopf, Phys. Rept. 858, 1 (2020), 1912.02253
  * (4) C. Egerer et al. (HadStruc), Phys. Rev. D 105, 034507 (2022), 2111.01808
  * (5) G. Backus, F. Gilbert, J. Roy. Astron. Sec 16, 196 (1968)
  * (6) M. Hansen, A. Lupo, N. Tantalo, Physical Review D 99, 094508 (2019)
  * (7) S.W. Mages, S. Borsányi, Z. Fodor, A. Schäfer, K. Szabó, _Shear viscosity from lattice QCD_ , in _The 32nd International Symposium on Lattice Field Theory_ (SISSA Medialab, 2015), Vol. 214, p. 232
  * (8) Y. Burnier, H.T. Ding, O. Kaczmarek, A.L. Kruse, M. Laine, H. Ohno, H. Sandmeyer, Journal of High Energy Physics 2017, 1 (2017)
  * (9) N.Y. Astrakhantsev, V.V. Braguta, A.Y. Kotov, Phys. Rev. D 98, 054515 (2018), 1804.02382
  * (10) A.N. Tikhonov, _On the stability of inverse problems_ , in _Dokl. Akad. Nauk SSSR_ (1943), Vol. 39, pp. 195–198
  * (11) M. Jarrell, J.E. Gubernatis, Physics Reports 269, 133 (1996)
  * (12) M. Asakawa, Y. Nakahara, T. Hatsuda, Progress in Particle and Nuclear Physics 46, 459 (2001)
  * (13) Y. Burnier, A. Rothkopf, Physical review letters 111, 182003 (2013)
  * (14) J. Skilling, in _Maximum-entropy and Bayesian methods in science and engineering_ (Springer, 1988), pp. 173–187
  * (15) G. Aarts, C. Allton, T. Harris, S. Kim, M.P. Lombardo, S.M. Ryan, J.I. Skullerud, JHEP 07, 097 (2014), 1402.6210
  * (16) S. Kim, P. Petreczky, A. Rothkopf, JHEP 11, 088 (2018), 1808.08781
  * (17) R. Larsen, S. Meinel, S. Mukherjee, P. Petreczky, Phys. Lett. B 800, 135119 (2020), 1910.07374
  * (18) Y. Burnier, O. Kaczmarek, A. Rothkopf, JHEP 12, 101 (2015), 1509.07366
  * (19) Y. Burnier, A. Rothkopf, Phys. Rev. D 95, 054511 (2017), 1607.04049
  * (20) D. Bala, O. Kaczmarek, R. Larsen, S. Mukherjee, G. Parkar, P. Petreczky, A. Rothkopf, J.H. Weber (HotQCD), Phys. Rev. D 105, 054513 (2022), 2110.11659
  * (21) J. Liang, T. Draper, K.F. Liu, A. Rothkopf, Y.B. Yang (XQCD), Phys. Rev. D 101, 114503 (2020), 1906.05312
  * (22) J. Karpie, K. Orginos, A. Rothkopf, S. Zafeiropoulos, JHEP 04, 057 (2019), 1901.05408
  * (23) S. Offler, G. Aarts, C. Allton, B. Jäger, S. Kim, M.P. Lombardo, B. Page, S.M. Ryan, J.I. Skullerud, T. Spriggs, PoS LATTICE2021, 509 (2022), 2112.02116
  * (24) L. Kades, J.M. Pawlowski, A. Rothkopf, M. Scherzer, J.M. Urban, S.J. Wetzel, N. Wink, F.P.G. Ziegler, Phys. Rev. D 102, 096001 (2020), 1905.04305
  * (25) R. Fournier, L. Wang, O.V. Yazyev, Q. Wu, Physical Review Letters 124, 056401 (2020)
  * (26) S. Shi, L. Wang, K. Zhou, Comput. Phys. Commun. 282, 108547 (2023), 2201.02564
  * (27) J. Horak, J.M. Pawlowski, J. Rodríguez-Quintero, J. Turnwald, J.M. Urban, N. Wink, S. Zafeiropoulos, Phys. Rev. D 105, 036014 (2022), 2107.13464
  * (28) R.A. Tripolt, P. Gubler, M. Ulybyshev, L. Von Smekal, Comput. Phys. Commun. 237, 129 (2019), 1801.10348
  * (29) C. Gattringer, K. Langfeld, International Journal of Modern Physics A 31, 1643007 (2016), number: 22 arXiv: 1603.09517
  * (30) J. Berges, K. Boguslavski, S. Schlichting, R. Venugopalan, Phys. Rev. D 89, 114007 (2014), 1311.3005
  * (31) G.D. Moore, M. Tassler, JHEP 02, 105 (2011), 1011.1167
  * (32) M. Cristoforetti, F. Di Renzo, L. Scorzato, AuroraScience Collaboration, Physical Review D 86, 074506 (2012), number: 7
  * (33) M. Namiki, I. Ohba, K. Okano, Y. Yamanaka, A.K. Kapoor, H. Nakazato, S. Tanaka, Lect. Notes Phys. Monogr. 9, 1 (1992)
  * (34) M. Laine, O. Philipsen, M. Tassler, JHEP 09, 066 (2007), 0707.2458
  * (35) A. Lehmann, A. Rothkopf, JHEP 07, 067 (2021), 2012.10089
  * (36) Y. Burnier, A. Rothkopf, Phys. Rev. D 86, 051503 (2012), 1208.1899
  * (37) K. Boguslavski, B.S. Kasmaei, M. Strickland, JHEP 10, 083 (2021), 2102.12587
  * (38) A. Rothkopf (2021), 2102.08616
  * (39) A. Rothkopf, J. Nordström (2022), 2205.14028
  * (40) D. Alvestad, R. Larsen, A. Rothkopf, JHEP 08, 138 (2021), 2105.02735
  * (41) G. Aarts, F.A. James, E. Seiler, I.O. Stamatescu, Physics Letters B 687, 154 (2010), adaptive step-size, 0912.0617
  * (42) D. Alvestad, R. Larsen, A. Rothkopf (2022)
  * (43) M. Troyer, U.J. Wiese, Phys. Rev. Lett. 94, 170201 (2005), cond-mat/0408370
  * (44) J. Berges, S. BorsÃ¡nyi, D. Sexty, I.O. Stamatescu, Physical Review D 75 (2007), hep-lat/0609058
  * (45) K. Boguslavski, P. Hotzy, D.I. Müller, _A stabilizing kernel for complex Langevin simulations of real-time gauge theories_ , in _39th International Symposium on Lattice Field Theory_ (2022), 2210.08020
  * (46) S. Woodward, P.M. Saffin, Z.G. Mou, A. Tranberg (2022), 2204.10101