11institutetext: Department of Physics, College of Science, Swansea
University, Swansea SA2 8PP, United Kingdom

# Electrical conductivity of the quark-gluon plasma: perspective from lattice
QCD

Gert Aarts Email<EMAIL_ADDRESS>Aleksandr Nikolaev Email:
<EMAIL_ADDRESS>
(August 27, 2020 )

###### Abstract

A discussion on the electrical conductivity of the quark-gluon plasma as
determined by lattice QCD is given. After a reminder of basic definitions and
expectations, various methods for spectral reconstruction are reviewed,
including the use of Ansätze and sum rules, the Maximum Entropy and Backus-
Gilbert methods, and Tikhonov regularisation. A comprehensive overview of
lattice QCD results obtained so far is given, including a comparison of the
different lattice formulations. A noticeable consistency for the
conductivities obtained is seen, in spite of the differences in the lattice
setups and spectral reconstruction methods. It is found that in the case of
quenched QCD little temperature dependence of $\sigma/T$ is seen in the
temperature range investigated, while for QCD with dynamical quarks a
reduction of $\sigma/T$ in the vicinity of the thermal crossover is observed,
compared to its value in the QGP. Several open questions are posed at the end.

Invited contribution to the EPJA topical issue “Theory of hot matter and
relativistic heavy-ion collisions (THOR)”.

###### pacs:

12.38.MhQuark-gluon plasma and 12.38.GcLattice QCD calculations

##  1 Introduction

The experimental heavy-ion programmes at the Relativistic Heavy Ion Collider
(RHIC) and the Large Hadron Collider (LHC), and, in the near future, at the
Nuclotron-based Ion Collider fAcility (NICA) and the Facility for Antiproton
and Ion Research (FAIR), offer exciting probes into the dynamics of strongly
interacting matter under extreme conditions. The relation with the underlying
theory, Quantum Chromodynamics, is established via phenomenology, which
permits a connection between quantities computable from first principles, such
as the equation of state, and measurable observables in the experiments.

In this contribution we focus on one such quantity, the electrical
conductivity $\sigma$ of the quark-gluon plasma (QGP). As we will review in
the next section, the usual definition of the conductivity, employing the Kubo
formula, relates it to a specific limit of Green’s functions in quantum field
theory, allowing for a computational formulation using, e.g., lattice QCD in
principle. On the other hand, the conductivity plays a role in charge
transport, particle production and the time evolution of electromagnetic
fields generated in heavy-ion collisions, see e.g. Refs. Deng:2012pc ;
Hirono:2012rt ; Tuchin:2013ie ; Yin:2013kya ; Kharzeev:2015znc ; Pratt:2019pnd
and references therein, emphasising its phenomenological relevance.

On the theoretical side, the conductivity can be computed using a variety of
methods, ranging from Feynman diagrams at weak coupling Aarts:2002tn ;
Gagnon:2006hi  and kinetic theory in QCD Arnold:2000dr ; Arnold:2003zc  or
effective models Cassing:2013iz ; Marty:2013ita  to holographic methods at
strong coupling Policastro:2002se ; Teaney:2006nc . Here we will focus on the
results obtained using numerical simulations of QCD discretised on the
lattice, as a first-principle tool to access nonperturbative information in
the vicinity of the deconfinement transition.

So far there are ${\cal O}(10)$ papers which have attempted to compute the
conductivity on the lattice Gupta:2003zh ; Aarts:2007wj ; Ding:2010ga ;
Francis:2011bt ; Brandt:2012jc ; Amato:2013naa ; Aarts:2014nba ;
Brandt:2015aqk ; Ding:2016hua ; Astrakhantsev:2019zkr . These papers differ
substantially in detail, partly indicating the increase in available computing
power over the past 15 years or so. For instance, there are simulations with
$N_{f}=0$ flavours (quenched QCD), and $N_{f}=2$ and $2+1$ dynamical flavours;
with quarks heavier than in nature or at the physical point; using a continuum
extrapolation or at fixed lattice spacing; with isotropic and anisotropic
($a_{\tau}\ll a_{s}$) lattices, etc (a detailed comparison is given in Sec.
4). Importantly, the methods used to extract the conductivity from the
Euclidean lattice correlators differ substantially as well, and include the
use of Ansätze, Bayesian approaches such as the Maximum Entropy Method, and
other regularisations. Despite these differences, a consistent picture is seen
to emerge, with approximate agreement between simulations with either
dynamical quarks or in the quenched case. The aim of this review is to give a
comprehensive overview of what has been obtained so far and provide a
comparison of the results. For completeness, we note here that we restrict
ourselves to the conductivity in the case of light quarks; we will not discuss
heavy-quark diffusion Meyer:2010tt ; Banerjee:2011ra ; Francis:2015daa  and
neither other transport coefficients such as the shear Meyer:2007ic ;
Astrakhantsev:2017nrs ; Pasztor:2018yae  and bulk viscosities Meyer:2007dy ;
Astrakhantsev:2018oue .

This paper is organised as follows. In the following section we present some
basic expressions relating the conductivity to various Green’s functions,
notably the spectral function and the corresponding Euclidean correlator,
using the Kubo relation. Some general remarks on expectations at high
temperature and the so-called transport peak are given as well. In Section 3
we discuss the various approaches that have been employed to reconstruct the
spectral function and extract the conductivity, given a numerically determined
Euclidean correlator. An overview of available lattice results is given in
Section 4, including a comparison between the values of $\sigma$ obtained so
far. Some related developments are summarised in Section 5. The final section
contains a summary, including some open questions. In the Appendix some well-
known relations between the various Green’s functions are collected.

##  2 Kubo formula and spectral function

The electromagnetic current in QCD receives contributions from all quark
flavours and reads

$j_{\mu}^{\rm em}(x)=\sum_{f=1}^{N_{f}}(eq_{f})j_{\mu}^{f}(x),\quad
j_{\mu}^{f}(x)=\bar{\psi}^{f}(x)\gamma_{\mu}\psi^{f}(x).$ (1)

Here $q_{f}$ denotes the fractional charge of the quark ($2/3$ or $-1/3$) and
$e$ the elementary charge. We restrict the discussion to light quarks, with
$N_{f}=2$ or $2+1$. The current is hermitian, ${j_{\mu}^{\rm
em}}^{\dagger}(x)=j_{\mu}^{\rm em}(x)$.

The electrical conductivity $\sigma$ indicates the linear relationship between
the current density and an electric field, $j_{i}^{\rm em}=\sigma E_{i}$,
according to Ohm’s law. Using linear-response theory, it can be related to the
current-current correlator in thermal equilibrium, in absence of the external
electric field, see e.g. Ref. Arnold:2000dr . More precisely, the conductivity
is proportional to the slope of the current-current spectral function in
thermal equilibrium,

$\rho_{\mu\nu}^{\rm em}(\omega,{\mathbf{p}})=\int d^{4}x\,e^{i\omega
t-{\mathbf{p}}\cdot{\mathbf{x}}}\langle[j_{\mu}^{\rm
em}(t,{\mathbf{x}}),j_{\nu}^{\rm em}(0,{\mathbf{0}})]\rangle,$ (2)

at vanishing energy and momentum, i.e.,

$\sigma=\frac{1}{6}\frac{\partial}{\partial\omega}\rho^{\rm
em}_{ii}(\omega,{\mathbf{0}})\Big{|}_{\omega=0}.$ (3)

Here the summation over spatial components, $i=1,2,3$, is understood. The
current-current spectral function is the expectation value of the commutator
of the electromagnetic current, evaluated at temperature $T$.

The conductivity is closely related Arnold:2000dr  to the charge diffusion
coefficient $D$, according to the Einstein relation,

$\sigma=\chi_{Q}D,$ (4)

where $\chi_{Q}$ is the charge susceptibility,

$\chi_{Q}=\frac{1}{TV}\left\langle(Q-\langle Q\rangle)^{2}\right\rangle.$ (5)

Here $V$ is the spatial volume and $Q$ is the total charge, i.e. the volume
integral of $j_{0}^{\rm em}(x)$.

Some well-known relations between the spectral function and other Green’s
functions are given in Appendix A. In particular, the spectral function is
related to the Euclidean correlator,

$G_{\mu\nu}^{\rm em}(\tau,{\mathbf{x}})=\langle j_{\mu}^{\rm
em}(\tau,{\mathbf{x}})j_{\nu}^{\rm em}(0,{\mathbf{0}})\rangle,$ (6)

via the standard relation [see Eq. (69)]

$G_{\mu\nu}^{\rm
em}(\tau,{\mathbf{p}})=\int_{0}^{\infty}\frac{d\omega}{2\pi}\,K(\tau,\omega)\rho_{\mu\nu}^{\rm
em}(\omega,{\mathbf{p}}),$ (7)

with the kernel

$K(\tau,\omega)=\frac{\cosh[\omega(\tau-1/2T)]}{\sinh(\omega/2T)}.$ (8)

Here the Euclidean time $0\leq\tau<1/T$. The question of computing the
conductivity on the lattice therefore boils down to numerically computing Eq.
(6), inverting Eq. (7) and extracting the slope according to Eq. (3). From now
we work at vanishing spatial momentum and drop the ${\mathbf{p}}$ dependence.
When prefactors involving $eq_{f}$ are dropped, the superscript ‘em’ is
omitted as well.

Figure 1: One-loop contribution to the current-current spectral function.

In order to prepare for the discussion of the lattice QCD results below, it is
useful to recall what can be expected at very high temperature, where QCD is
weakly coupled. At leading (zeroth) order in perturbation theory (see Fig. 1),
the spectral function, for a single flavour with mass $m$, reads Aarts:2005hg
111Note there is a typo in the last line of Eq. (19) in Ref. Aarts:2005hg :
$a^{(1)}_{H}+a^{(2)}_{H}$ should read $a^{(1)}_{H}+a^{(3)}_{H}$.

$\displaystyle\rho_{ii}(\omega)$ $\displaystyle=$ $\displaystyle 2\pi
N_{c}I\omega\delta(\omega)$ (9)
$\displaystyle+\frac{N_{c}}{2\pi}\theta(\omega^{2}-4m^{2})\sqrt{\frac{\omega^{2}-4m^{2}}{\omega^{2}}}$
$\displaystyle\times\left(\omega^{2}+2m^{2}\right)\left[1-2n_{F}\left(\frac{\omega}{2}\right)\right],$

with $N_{c}=3$. Here $n_{F}(\omega)=1/[\exp(\omega/T)+1]$ is the Fermi-Dirac
distribution. The quantity $I$ in the first term reads

$I=-4\int\frac{d^{3}k}{(2\pi)^{3}}\,n_{F}^{\prime}(\omega_{\mathbf{k}})\frac{k^{2}}{\omega_{\mathbf{k}}^{2}},$
(10)

with $\omega_{\mathbf{k}}=\sqrt{{\mathbf{k}}^{2}+m^{2}}$. For massless quarks
this evaluates as

$I\big{|}_{m=0}=\frac{T^{2}}{3}.$ (11)

Since

$1-2n_{F}\left(\frac{\omega}{2}\right)=\tanh\left(\frac{\omega}{4T}\right),$
(12)

the spectral function is odd, $\rho_{ii}(-\omega)=-\rho_{ii}(\omega)$, as it
should be. Below we take $\omega\geq 0$. The corresponding Euclidean
correlator reads, in the massless limit Aarts:2005hg ,

$G_{ii}(\tau)=N_{c}T^{3}\left[\frac{1}{3}+\frac{3u+u\cos(2u)-2\sin(2u)}{\sin^{3}(u)}\right],$
(13)

where $u=2\pi T(\tau-1/2T)$. The first (constant) term comes from the first
term in Eq. (9); the $\tau$ dependent term from the second one. At the
midpoint, $\tau=1/2T$, the contributions from both terms are comparable,
$G_{ii}(\tau=1/2T)=N_{c}T^{3}(1/3+2/3)=N_{c}T^{3}$.

Let us now discuss the two contributions in Eq. (9) in more detail. We start
with the second one. This term arises from the cut of the one-loop
polarisation diagram in Fig. 1 corresponding to a decay process, with
$\omega=\omega_{\mathbf{k}}+\omega_{{\mathbf{p}}+{\mathbf{k}}}$ (with
${\mathbf{k}}$ the loop momentum and ${\mathbf{p}}={\mathbf{0}}$), which is
permissible once $|\omega|>|2m|$. It contains the vacuum contribution,
increasing as $\omega^{2}$ at large $\omega$, and a thermal contribution,
leading to Pauli blocking, which is exponentially suppressed at large
energies. Note that this term does not contribute to the conductivity; in the
massless limit it increases as $\omega^{3}$ as $\omega\to 0$. We will refer to
this term as the continuum or perturbative contribution.

The first term, with $\omega\delta(\omega)$, is only present at nonzero
temperature and arises from the cut corresponding to scattering with particles
in the heatbath,
$\omega_{\mathbf{k}}+\omega=\omega_{{\mathbf{p}}+{\mathbf{k}}}$, in the limit
that ${\mathbf{p}}\to{\mathbf{0}}$. Extracting the conductivity from this term
yields infinity, reflecting the fact that for free particles the mean free
path and hence the conductivity diverges. Interactions make the mean free path
finite, due to scattering in the plasma. The result is that the $\delta$
function in Eq. (9) is smeared out and takes the form of a so-called transport
peak,

$2\pi N_{c}I\omega\delta(\omega)\to\rho_{\rm trans}(\omega)=A_{\rm
trans}\frac{\gamma\omega}{\omega^{2}+\gamma^{2}},$ (14)

where $\gamma$ is proportional to the collisional scattering width or the
inverse mean free path, $A_{\rm trans}(=2N_{c}I)$ is the overall coefficient,
and $\sigma\propto A_{\rm trans}/\gamma$. The origin of the transport peak can
be seen in various ways, using e.g. Feynman diagrams and pinching poles
Aarts:2002cc ; Aarts:2002vx  or kinetic theory Petreczky:2005nh ; Meyer:2011gj
. It should be noted that this form of transport peak is the simplest form
encountered. An important observation is that in the case of a narrow
transport peak ($\gamma\ll T$), as is the case for weakly-coupled theories,
the Euclidean correlator is not sensitive to details of the transport peak but
only to its area Aarts:2002cc . In the limit that $\omega<\Lambda\ll T$, the
kernel (8) simplifies to $2T/\omega$, and integrating the transport peak (14)
in Eq. (7) yields

$\displaystyle G_{\rm
trans}(\tau)=\int_{0}^{\Lambda}\frac{d\omega}{2\pi}\,K(\tau,\omega)\rho_{\rm
trans}(\omega)$ $\displaystyle\sim A_{\rm
trans}\int_{0}^{\Lambda}\frac{d\omega}{2\pi}\,\frac{2T}{\omega}\frac{\gamma\omega}{\omega^{2}+\gamma^{2}}\sim\frac{1}{2}A_{\rm
trans}T,$ (15)

where in the last expression the cutoff $\Lambda$ on $\omega$ has been
removed. The crucial observation is that this expression is independent of
$\gamma$ and the Euclidean time $\tau$, indicating the insensitivity of the
correlator to narrow transport peaks. In fact, taking $A_{\rm trans}=2N_{c}I$,
its value is the same as in the non-interacting theory.

So far we have only considered the diagram given in Fig. 1, dressed with
gluons and closed loops of sea quarks in the presence of interactions (in
lattice QCD this is referred to as the connected contribution). However, when
interactions are included, there are also contributions from diagrams with a
different topology, namely with two closed fermion loops, connected via
gluons. These arise from Wick contractions of the current operator,
$j_{i}=\bar{\psi}\gamma_{i}\psi$, with itself. In lattice QCD, this is
commonly referred to as the disconnected contribution. Perturbatively, they
start contributing at ${\cal O}(\alpha_{s}^{3})$ only, and hence are
suppressed at very high $T$. Another distinction between the connected and
disconnected contributions concerns the appearance of the electromagnetic
charges. Take for simplicity $N_{f}=2$ or $3$ degenerate flavours. In the
connected contribution one finds the sum over charges squared, which is
usually denoted as

$C_{\rm em}=\sum_{f}(eq_{f})^{2}=\begin{cases}\frac{5}{9}e^{2}\qquad(u,d),\\\
\frac{6}{9}e^{2}\qquad(u,d,s).\end{cases}$ (16)

In the disconnected contribution on the other hand, we find the square of the
sums,

$C_{\rm em}^{\rm
disc}=\big{(}\sum_{f}eq_{f}\big{)}^{2}=\begin{cases}\frac{1}{9}e^{2}\qquad(u,d),\\\
0\qquad\quad(u,d,s).\end{cases}$ (17)

Hence it is noted that the disconnected contribution vanishes for three
degenerate flavours.

Up to now we have focussed on the expectation at high temperature. At low
temperature, in the hadronic phase, the current-current correlator couples to
vector mesons, with details depending on the flavour content. A comprehensive
discussion in the $N_{f}=2$ case can be found in Ref. Brandt:2012jc . In the
spectral function one therefore expects bound-state peaks, representing
mesonic ground and excited states. As the system is heated, these bound states
are expected to dissolve and the high-temperature behaviour to emerge.

To conclude this section, we note that lattice QCD simulations include quarks
scattering with gluons, i.e. the electromagnetic field and other charge
carriers (leptons) are not included. Hence the conductivity reviewed here
indicates the contribution from the strong interaction only, yielding insight
into the strongly coupled nature of the quark-gluon plasma.

##  3 Spectral reconstruction

The main problem in extracting the conductivity from numerically determined
Euclidean lattice correlators is the inversion problem, see Eq. (7). As a
reminder, on the lattice temperature is encoded in the compact direction in
imaginary time, with circumference $1/T=a_{\tau}N_{\tau}$. Here $a_{\tau}$ is
the temporal lattice spacing and $N_{\tau}$ the number of points in the time
direction; Euclidean time is discretised as
$\tau/a_{\tau}=0,\ldots,N_{\tau}-1$. Since the correlator $G(\tau)$ is known
at a finite number of temporal points and the spectral function $\rho(\omega)$
is in principle a continuous function of $\omega$, this inversion problem is
far from straightforward. To be more precise, due to reflection symmetry,
$K(\tau,\omega)=K(1/T-\tau,\omega)$, the number of points available for the
analysis is on the order of $N_{\tau}/2$; even after placing an upper limit on
the $\omega$ interval, such that $0<\omega<\omega_{\rm max}$, and discretising
the finite interval, typically on the order of $N_{\omega}=1000$ points are
used to present $\rho(\omega)$. Since $N_{\omega}\gg N_{\tau}$, the inversion
problem is ill-posed. In addition, the focus on the $\omega\to 0$ limit makes
the inversion more challenging than for spectral functions in general, when
the interest is in frequencies on the order of the temperature or above, as
the discussion around Eq. (15) indicates.

Several methods have been developed to tackle this problem. Here we briefly
review the ones applied to the conductivity. It is fair to state that no
single method is yet fully robust on its own. Hence it is of interest to
compare and contrast the results obtained so far, and seek for
(in)consistencies. This will be done in Sec. 4.

###  3.1 Reconstructed correlators

Before investigating the temperature dependence of the spectral function, we
note that the Euclidean correlator (7) depends on temperature in two ways:

  * •

via the temperature dependence of the kernel, $K(\tau,\omega)$, due to the
compact time direction, $0\leq\tau<1/T$;

  * •

via the temperature dependence of the spectral function, due to changes in the
quark-gluon plasma.

The first effect leads to temperature dependence of the correlator even when
the spectral function is unchanged. It is important to disentangle this from
the sought second (physical) effect due to actual changes in the plasma. This
can be investigated using so-called reconstructed correlators Datta:2003ww ;
Aarts:2007pk . Let us suppose a spectral function $\rho(\omega;T_{0})$ is
determined (with some confidence) at a reference temperature $T_{0}$. Assuming
that the spectral function is unchanged, a correlator at a different
temperature $T$ can then be defined as

$G_{\rm
recon}(\tau,T;T_{0})=\int_{0}^{\infty}\frac{d\omega}{2\pi}\,K(\tau,\omega;T)\rho(\omega;T_{0}).$
(18)

This construction takes into account the trivial temperature dependence due to
the kernel, the first effect above. Comparing this reconstructed correlator
with the actual correlator at temperature $T$ allows one to draw conclusions
on the second effect, i.e. changes in the spectral function due to a change in
the physical situation. A difference between the actual and the reconstructed
correlator implies a change in the spectral function (the inverse is not
necessarily true).

###  3.2 Sum rules

Exact sum rules are important Brandt:2012jc  to constrain the current-current
spectral function at nonzero temperature. Defining the difference between the
finite and zero-temperature spectral function as

$\Delta\rho(\omega,T)\equiv\rho_{ii}(\omega,T)-\rho_{ii}(\omega,0),$ (19)

one finds the sum rule Brandt:2012jc , in the thermodynamic limit,

$\int_{-\infty}^{\infty}\frac{d\omega}{\omega}\Delta\rho(\omega,T)=0.$ (20)

Note that the zero-temperature $\omega^{2}$ contribution cancels in the
subtraction (19). Since the Operator Product Expansion (OPE) predicts
CaronHuot:2009ns  that the thermal contribution decays as $(T/\omega)^{2}$ at
large $\omega$, the integral in the sum rule converges. One may verify that
this sum rule indeed holds for free fermions, using Eq. (9).

The sum rule indicates that enhancement of spectral weight at small energies,
i.e. due to a larger transport peak, should be compensated by a loss of
spectral weight elsewhere. Since the sum rule is exact, it should be satisfied
by reconstructed spectral functions on the lattice, where it can be
implemented as a check or a constraint. This sum rule, and two additional
ones, are further analysed in Refs. Gubler:2016hnf ; Gubler:2017qbs .

###  3.3 Ansätze

The first step to resolve the ill-posedness of the inversion is to reduce the
number of parameters needed to model the spectral function. The easiest way to
do so is by providing an Ansatz for $\rho(\omega)$ with less fit parameters
than data points. The downside is that this introduces an obvious bias, which
is difficult to avoid. Moreover, since the spectral function is expected to
behave in quite a different manner in the low- and high-temperature phases,
the Ansatz has to be sufficiently rich to capture this. Some features to be
included are

  * •

a transport peak at small $\omega$, with in particular a linear slope in
$\omega$;

  * •

continuum ($\omega^{2}$) contribution at high $\omega$, possibly modified by
lattice artefacts Karsch:2003wy ; Aarts:2005hg ;

  * •

at least one bound-state peak in the low-temperature phase, to represent the
vector meson.

Refs. Ding:2010ga ; Francis:2011bt ; Ding:2016hua  employ an Ansatz combining
a transport peak and the expected perturbative continuum behaviour (for
massless quarks) in the deconfined phase,

$\rho(\omega)=\rho_{\rm trans}(\omega)+\rho_{\rm pert}(\omega),$ (21)

where $\rho_{\rm trans}(\omega)$ is given in Eq. (14) and

$\rho_{\rm pert}(\omega)=\frac{3}{2\pi}A_{\rm
pert}\omega^{2}\left[1-2n_{F}\left(\frac{\omega}{2}\right)\right],$ (22)

c.f. Eq. (9). The three temperature-dependent parameters are the coefficients
$A_{\rm trans,pert}$ and the width $\gamma$ of the transport peak. Note that
$A_{\rm pert}=1$ for free fermions; it parametrizes deviations from a free
spectral function at large energies. As stated, the functional form of this
Ansatz is the combination of two functions. Modifying the transport peak to a
flat featureless function, as seen e.g. in holography Teaney:2006nc , Ref.
Ding:2016hua  finds that the data may not have the resolution to differentiate
between these two shapes. This is incorporated in the systematic uncertainty
of the final quoted result for $\sigma$ Ding:2016hua .

Ref. Brandt:2012jc  employs a related Ansatz for the subtracted spectral
function, $\Delta\rho(\omega,T)$, defined in Eq. (19). In addition to the
transport peak and the continuum contribution, this Ansatz also includes a
bound-state peak. It reads

$\Delta\rho(\omega)=\rho_{\rm trans}(\omega)+\Delta\rho_{\rm
pert}(\omega)+\rho_{\rm bound}(\omega),$ (23)

with the new term

$\rho_{\rm bound}(\omega)=A_{\rm
bound}\frac{2g_{B}\tanh(\omega/T)^{3}}{4(\omega-m_{B})^{2}+g_{B}^{2}}.$ (24)

Here $m_{B},g_{B}$ and $A_{\rm bound}$ indicate the mass, width and strength
of the bound state. The factor $\tanh(\omega/T)^{3}$ ensures the contribution
does not contribute to the conductivity in the $\omega\to 0$ limit and decays
as $1/\omega^{2}$ at large $\omega$, as predicted by the OPE CaronHuot:2009ns
. It is also noted in Refs. Burnier:2012ts ; Brandt:2012jc  that the transport
peak (14) in fact violates this condition. Hence it is proposed Brandt:2012jc
to modify it as

$\rho_{\rm trans,mod}(\omega)=A_{\rm
trans}\frac{T\tanh(\omega/T)\gamma}{\omega^{2}+\gamma^{2}}\,,$ (25)

which still has linear behaviour at small $\omega$ but decays as
$1/\omega^{2}$ at large $\omega$. Finally, the main reason for the subtraction
is to eliminate the zero-temperature $\omega^{2}$ contribution,

$\Delta\rho_{\rm pert}(\omega)=\rho_{\rm pert}(\omega;T)-\rho_{\rm
pert}(\omega;T=0),$ (26)

which eliminates the “1” in Eq. (22). This allows the analysis to focus on
frequencies on the order of the temperature, without being overwhelmed by the
$\omega^{2}$ term. Overall, the number of parameters ($A_{\rm
trans,pert,bound},m_{B},g_{B},\gamma$) to be fitted is quite large, which is
carried out by fixing them in steps, while satisfying the sum rule (20). A
reduced model, using only $\rho_{\rm trans}+\rho_{\rm pert}$, is employed as
well. A variation of this Ansatz is also used in Ref. Brandt:2015aqk ,
replacing the bound state by a delta-function and introducing explicit
thresholds for the various terms.

###  3.4 Maximum Entropy Method

The Ansätze described above have to incorporate a wide range of physics input
(bound states, transport peak, continuum contribution), each of which is
defined by a number of parameters, making the fit highly nonlinear and
depending on the choice of model functions. It is therefore desirable to use
model-independent reconstruction methods for the spectral function. It will be
necessary to regularise standard minimisation procedures, due to the ill-
posedness of the inversion. Before proceeding, we note the possibility to
rescale the kernel and the spectral function,

$K(\tau,\omega)\to
f(\omega)K(\tau,\omega),\qquad\rho(\omega)\to\rho(\omega)/f(\omega),$ (27)

leaving the product unchanged, to stabilise the inversion. A common rescaling
is to use $f(\omega)=\omega$, to resolve the $1/\omega$ divergence in the
kernel as $\omega\to 0$ Aarts:2007wj .

There is still a requirement to reduce the number of parameters to be
determined, to make the inversion solvable. There is considerable freedom to
do so. For instance, in the Maximum Entropy Method (MEM), and in particular
Bryan’s method Bryan , the (rescaled) spectral function is parametrised as
Aarts:2007wj ; Asakawa:2000tr

$\frac{\rho(\omega)}{\omega}=\frac{m(\omega)}{\omega}\exp\sum_{i=1}^{N_{\rm
coeff}}c_{i}u_{i}(\omega),$ (28)

with $u_{i}(\omega)$ $(i=1,\ldots,N_{\rm coeff})$ an orthogonal but incomplete
set of basis functions,

$\int_{0}^{\infty}\frac{d\omega}{2\pi}\,u_{i}(\omega)u_{j}(\omega)=\delta_{ij}.$
(29)

The reduction follows since $N_{\rm coeff}\sim N_{\tau}/2\ll N_{\omega}$. The
form (28) is motivated by positivity, $\rho(\omega)/\omega\geq 0$. In MEM,
$m(\omega)$ is referred to as the default model, see below. The conductivity
is now determined by

$\sigma\sim m^{\prime}(0)\exp\sum_{i}c_{i}u_{i}(0),$ (30)

and hence depends on all coefficients and the default model.

In MEM the coefficients $c_{i}$ are determined by constructing the most
probable spectral function as defined by the extremum of the conditional
probability $P(\rho|DH)$ Asakawa:2000tr . Here $D$ indicates the data and $H$
additional prior knowledge. The method relies on Bayes’ theorem,

$P(\rho|DH)=\frac{P(D|\rho H)P(\rho|H)}{P(D|H)},$ (31)

where $P(A|B)$ stands for the conditional probability of $A$ given $B$. In
this expression, $P(D|\rho H)$ is the likelihood function, $P(\rho|H)$ the
prior probability, and $P(D|H)$ a normalisation. While the likelihood
function, $P(D|\rho H)=e^{-L(\rho)}$, is familiar from standard
$\chi^{2}$-minimisation, the prior probability contains an entropy-like term,

$P(\rho|H)=e^{\alpha S(\rho)},$ (32)

with

$S(\rho)=\int_{0}^{\infty}\frac{d\omega}{2\pi}\,\left[\rho(\omega)-m(\omega)-\rho(\omega)\ln\frac{\rho(\omega)}{m(\omega)}\right],$
(33)

giving the method its name. The conditional probability now reads

$P(\rho|DH)\propto e^{-L(\rho)+\alpha S(\rho)},$ (34)

with $\alpha$ determining the balance between the two terms. While at
$\alpha=0$ the method reduces to a standard fitting procedure, in absence of
any data the probability is extremised when $\rho(\omega)=m(\omega)$, yielding
the default result. For further details on MEM we refer to Ref. Asakawa:2000tr
.

MEM has been applied to the conductivity in Refs. Gupta:2003zh ; Aarts:2007wj
; Ding:2010ga ; Francis:2011bt ; Amato:2013naa ; Aarts:2014nba . We note here
that the basis functions $u_{i}(\omega)$ are obtained via a singular value
decomposition (SVD) of the kernel $K(\tau,\omega)$, when viewed as a $N_{\rm
coeff}\times N_{\omega}$ matrix, where $N_{\rm coeff}\lesssim N_{\tau}/2$,
linking the size of the set of basis functions to the temporal extent. This is
a limitation for small $N_{\tau}$ and has motivated the use of anisotropic
lattices to increase the number of data points and basis functions at a given
temperature, see Ref. Amato:2013naa  and especially Ref. Aarts:2014nba  for a
systematic study. While at first sight the choice of default model appears to
play an important role in the formulation, in practice it has been found that
the dependence on $m(\omega)$ is quite mild Aarts:2014nba . The main
systematic uncertainty enters via the formulation of the method, such as the
choice of the prior probability, which can be addressed e.g. by a comparison
with other independent approaches.

###  3.5 Backus-Gilbert method

The Backus-Gilbert method is such an independent approach, designed for
solving linear ill-defined problems with controllable regularisation and
systematic uncertainty. Rather than reconstructing the entire spectral
function $\rho(\omega)$, it aims to represent it by an estimator

$\displaystyle\widehat{\rho}(\omega_{0})=\int_{0}^{\infty}d\omega\,\delta(\omega_{0},\omega)\rho(\omega),$
(35)

where $\delta(\omega_{0},\omega)$ is called the resolution function, which
should be narrowly peaked around $\omega_{0}$ and normalised,

$\int_{0}^{\infty}d\omega\,\delta(\omega_{0},\omega)=1,$ (36)

similar to a delta function. Ideally one wants to make the resolution function
as narrow as possible, for given correlator and kernel. For this purpose the
following linear Ansatz is assumed Brandt:2015aqk ; Astrakhantsev:2019zkr

$\delta(\omega_{0},\omega)=\sum_{n=1}^{N_{\tau}-1}q_{n}(\omega_{0})K(\tau_{n},\omega).$
(37)

The functions $q_{n}(\omega_{0})$ are found by minimising the second moment of
the resolution function squared,

$\Gamma_{\omega_{0}}=\int_{0}^{\infty}d\omega\,(\omega-\omega_{0})^{2}\delta^{2}(\omega_{0},\omega),$
(38)

which should effectively minimise its width. Combining the equations above,
one obtains the solution (summation over repeated indices implied)

$q_{n}(\omega_{0})=\frac{W(\omega_{0})^{-1}_{nm}R_{m}}{R_{k}W(\omega_{0})^{-1}_{kl}R_{l}},$
(39)

in terms of

$\displaystyle W(\omega_{0})_{nm}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega\,(\omega-\omega_{0})^{2}K(\tau_{n},\omega)K(\tau_{m},\omega),$
(40) $\displaystyle R_{n}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega\,K(\tau_{n},\omega)\,.$ (41)

Note that since the kernel $K(\tau,\omega)$ is symmetric around the midpoint,
the number of independent basis functions $q_{n}$ is limited by $N_{\tau}/2$.
In general, the larger the number of available time slices, the narrower
resolution functions can be constructed.

The spectral function can now be estimated by combining Eqs. (35, 37) and the
definition of the correlator (7), as

$\displaystyle\widehat{\rho}(\omega_{0})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega\,\sum_{n}q_{n}(\omega_{0})K(\tau_{n},\omega)\rho(\omega)$
(42) $\displaystyle=$ $\displaystyle
2\pi\sum_{n}q_{n}(\omega_{0})G(\tau_{n}).$

The conductivity is then extracted in a usual way,

$\sigma\sim\lim_{\omega_{0}\rightarrow
0}\frac{\widehat{\rho}(\omega_{0})}{\omega_{0}},$ (43)

assuming that the estimator $\widehat{\rho}(\omega_{0}\sim 0)$ is close to the
physical spectral function. The goal is therefore to find the optimal set of
functions $q_{n}(\omega_{0})$ which make $\delta(\omega_{0},\omega)$ as narrow
as possible in $\omega$ for $\omega_{0}\sim 0$.

The main difficulty in this formulation is the inversion of the matrix
$W(\omega_{0})$, which is usually ill-conditioned, due to the exponential
decay of the kernel. This can be ameliorated by rescaling, see Eq. (27),
changing the estimator to

$\widehat{\rho}(\omega_{0})=f(\omega_{0})\int_{0}^{\infty}d\omega\,\delta(\omega_{0},\omega)\frac{\rho(\omega)}{f(\omega)},$
(44)

and finding an optimal choice for $f(\omega)$ (e.g. $f(\omega)=\omega$).
Secondly, the matrix $W$ can be regularised as

$W_{nm}\rightarrow\lambda W_{nm}+(1-\lambda)S_{nm},$ (45)

where $S_{nm}$ is the covariance matrix for the correlator $G(\tau_{n})$, and
$\lambda$ is a tunable parameter, determined by comparing the behaviour of
$\delta(\omega_{0},\omega)$ for different values of $\lambda$. Further
discussion on successes and limitations of this approach in the context of the
conductivity can be found in Refs. Brandt:2015aqk ; Astrakhantsev:2019zkr .

###  3.6 Tikhonov regularisation

The method of additive regularisation, see Eq. (45), is not the only one.
Tikhonov regularisation ATikhonov:reg  acts as a complimentary instrument to
other reconstruction methods, where an inversion of an ill-conditioned matrix
has to be performed. Let us consider the $N\times N$ matrix $W$, defined in
Eq. (40), for which a straightforward inversion fails. A singular value
decomposition yields

$\displaystyle W$ $\displaystyle=$ $\displaystyle U\Sigma V^{T},\quad
U^{T}U=V^{T}V=1\\!\\!1,$ (46) $\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\mbox{diag}(\sigma_{1},\,\sigma_{2},\,\ldots,\,\sigma_{N}),$

with $\sigma_{1}\geq\sigma_{2}\geq\ldots\geq\sigma_{N}$. Since
$W^{-1}=V\Sigma^{-1}U^{T}$, the inversion of $W$ comes down to the inversion
of $\Sigma$, which can be regularised as follows,

$\Sigma^{-1}=\mbox{diag}\left(\frac{\sigma_{1}}{\sigma_{1}^{2}+\epsilon^{2}},\,\frac{\sigma_{2}}{\sigma_{2}^{2}+\epsilon^{2}},\,\ldots,\,\frac{\sigma_{N}}{\sigma_{N}^{2}+\epsilon^{2}}\right),$
(47)

where the parameter $\epsilon$ has to be chosen carefully; small $\epsilon$’s
lead to precise but unstable results, while large $\epsilon$’s guarantee
stable inversion at a cost of loss of accuracy. Further experiments with this
approach can be found in Ref. Astrakhantsev:2019zkr .

###  3.7 Other approaches

Here we briefly list some additional inversion methods. Refs. Burnier:2011jq ;
Burnier:2012ts  further develop the proposal of Ref. Cuniberti , which
formulates a unique analytic continuation which can be constructed explicitly,
provided the correlator satisfies certain asymptotic behaviour in Minkowski
time. The crucial requirement is that a continuum-extrapolated result for the
lattice correlator is available, and that short-distance divergences, present
at zero temperature, have been subtracted. So far these requirements are only
met in quenched QCD; in Sec. 4 the application of this approach will be
discussed further.

The following methods have not been yet been applied to the determination of
the conductivity, or other transport coefficients, in QCD, as far as we know.
The Maximum Entropy Method is only one of a number of Bayesian methods;
alternative Bayesian approaches can be found in Refs. Burnier:2013nla ;
Rothkopf:2011ef . Ref. Tripolt:2018xeo  proposes an inversion based on the
Schlessinger point or Resonances Via Padé method, which is based on a
rational-fraction representation similar to Padé approximation methods. It is
found that the method is competitive to MEM and Backus-Gilbert, provided the
errors of the input data are small enough. Interestingly, Ref. Tripolt:2018xeo
applies the method to the extraction of the conductivity in graphene,
described by a tight-binding model.

A very recent development is the implementation of machine learning approaches
to tackle spectral reconstruction. Supervised learning of fully connected as
well as convolutional neural networks was applied to mock data in Ref.
LKades2019spectral . In Ref. Arsenault  kernel ridge regression was applied to
mock data in quantum many-body physics; a first application to QCD data can be
found in Ref. Offler:2019eij  for bottomonium correlators. More developments
in the realm of machine learning are expected in the near future.

##  4 Lattice QCD results

In this section we give an overview of results obtained for the electrical
conductivity in QCD, with gauge group SU(3), for $N_{f}=0$ (quenched QCD) and
$N_{f}=2$ and $2+1$ dynamical flavours. We do not discuss results obtained in
effective models or in other gauge theories; see Sec. 5 for some results in
the SU(2) theory. The papers we discuss are listed in Table 1 in chronological
order, with some details on the ensembles. Refs. Gupta:2003zh ; Aarts:2007wj ;
Ding:2010ga ; Francis:2011bt ; Ding:2016hua  concern quenched QCD, while in
Refs. Brandt:2012jc ; Amato:2013naa ; Aarts:2014nba ; Brandt:2015aqk ;
Astrakhantsev:2019zkr  dynamical quarks are included ($N_{f}=2$ and $2+1$). In
the dynamical studies the sea quarks are of the Wilson-clover type, with a
pion mass heavier than in nature. The exception is Ref. Astrakhantsev:2019zkr
, with staggered sea quarks and a physical pion mass. All studies employ an
isotropic lattice, with $a_{s}=a_{\tau}$ (the spatial and temporal lattice
spacing respectively), except Refs. Amato:2013naa ; Aarts:2014nba , in which
anisotropic lattices, with $a_{s}/a_{\tau}=3.5$, are employed. The advantage
of the latter is the finer temperature resolution, when $N_{\tau}$ is varied.
The lattice cutoff is dealt with in a variety of ways:

  * •

continuum limit, indicated by $a_{\tau}\to 0$;

  * •

fixed scale: the lattice cutoff is fixed and temperature is varied by changing
$N_{\tau}$, according to $T=1/(a_{\tau}N_{\tau})$;

  * •

fixed cutoff: one lattice spacing is available at each temperature, lattice
spacings at different temperatures are not identical.

In simulations with dynamical quarks, continuum limits are not yet available
and the (temporal) lattice spacings lie in the range $0.0350<a_{\tau}<0.0618$
fm.

Table 1: Details of the lattice QCD ensembles to compute the electrical conductivity. Here $a_{\tau}$ and $a_{s}$ denote the temporal and spatial lattice spacing respectively. ref. | arXiv number |  $N_{f}$ (sea) | fermion type |  $m_{\pi}$ [MeV] |  $a_{\tau}$[fm] | $a_{s}/a_{\tau}$ | discretisation
---|---|---|---|---|---|---|---
Gupta:2003zh  | hep-lat/0301006 | 0 | quenched | $-$ | $a_{\tau}\to 0$ | 1 | continuum limit
Aarts:2007wj  | hep-lat/0703008 | 0 | quenched | $-$ | $0.0488,0.0203$ | 1 | fixed cutoff
Ding:2010ga  | 1012.4963 | 0 | quenched | $-$ | $a_{\tau}\to 0$ | 1 | continuum limit
Francis:2011bt  | 1112.4802 | 0 | quenched | $-$ | $0.015$ | 1 | fixed scale
Brandt:2012jc  | 1212.4200 | 2 | Wilson-clover | 270 | $0.0486(4)(5)$ | 1 | fixed cutoff
Amato:2013naa  | 1307.6763 | $2+1$ | Wilson-clover | 384(4) | $0.0350(2)$ | 3.5 | fixed scale
Aarts:2014nba  | 1412.6411 | $2+1$ | Wilson-clover | 384(4) | $0.0350(2)$ | 3.5 | fixed scale
Brandt:2015aqk  | 1512.07249 | $2$ | Wilson-clover | 270 | $0.0486(4)(5)$ | 1 | fixed scale
Ding:2016hua  | 1604.06712 | 0 | quenched | $-$ | $a_{\tau}\to 0$ | 1 | continuum limit
Astrakhantsev:2019zkr  | 1910.08516 | $2+1$ | staggered | 134.2(6) | $0.0618,0.0493$ | 1 | fixed cutoff
Table 2: Lattice sizes used compute the electrical conductivity: number of spatial ($N_{s}$) and temporal ($N_{\tau}$) lattice points, and corresponding temperatures. Details of the so-called “zero-temperature” lattices used for tuning are not listed. ref. | $N_{s}$ | $N_{\tau}$ | temperature
---|---|---|---
Gupta:2003zh  | 18, $\ldots$, 44 | 14, 12, 10, 8 | $T/T_{c}=1.5,2,3$
Aarts:2007wj  | 48, 64 | 24, 16 | $T/T_{c}=0.62,1.5,2.25$
Ding:2010ga  | 128 | 48, 32, 24, 16 | $T/T_{c}=1.45$
Francis:2011bt  | 128 | 40, 32, 16 | $T/T_{c}=1.16,1.49,2.98$
Brandt:2012jc  | 64 | 16 |  $T=250$ MeV
Amato:2013naa ; Aarts:2014nba  | 24, 32 | 48, 40, 36, 32, 28, 24, 20, 16 |  $T=117,141,156,176,201,235,281,352$ MeV
Brandt:2015aqk  | 64 | 24, 20, 16, 12 |  $T=169,203,254,338$ MeV
Ding:2016hua  | 96, 128, 144, 192 | 64, 56, 48, 42, 32, 28, 24 | $T/T_{c}=1.1,1.3,1.5$
Astrakhantsev:2019zkr  | 48, 64 | 10, 16 |  $T=200,250$ MeV
Table 3: Details of the current and inversion method used to compute the electrical conductivity. ref. | fermion type | current | renormalised | inversion method
---|---|---|---|---
Gupta:2003zh  | staggered | local | $-$ | Bayesian priors, MEM, Ansatz
Aarts:2007wj  | staggered | local | $-$ | MEM
Ding:2010ga  | Wilson-clover | local | ✓ | Ansatz, MEM
Francis:2011bt  | Wilson-clover | local | ✓ | Ansatz, MEM
Brandt:2012jc  | Wilson-clover | local | ✓ | Ansatz, sum rule constraints
Amato:2013naa ; Aarts:2014nba  | Wilson-clover | conserved | ✓ | MEM
Brandt:2015aqk  | Wilson-clover | mixed local-conserved | ✓ | Ansatz, sum rule constraints, Backus-Gilbert
Ding:2016hua  | Wilson-clover | local | ✓ | Ansatz
Astrakhantsev:2019zkr  | staggered | conserved | ✓ | Backus-Gilbert, Tikhonov regularisation

Table 2 contains the details of the lattice geometry, i.e. the number of
lattice points in spatial ($N_{s}$) and temporal ($N_{\tau}$) direction. The
corresponding temperatures are expressed in units of $T_{c}$ for the quenched
ensembles and in units of MeV for the dynamical ones. The largest number of
temperatures is considered in Refs. Amato:2013naa ; Aarts:2014nba , with 8
temperature values, ranging from 117 to 352 MeV.

Table 3 finally lists some details on the currents and methods used to compute
the conductivity:

  * •

so far the type of valence quarks used in the current have been either
staggered or Wilson-clover fermions. For staggered quarks, the current-current
correlator (6) has an oscillating structure, of the form Aarts:2007wj ;
Astrakhantsev:2019zkr ,

$\displaystyle G_{ij}(x)$ $\displaystyle=$ $\displaystyle\langle
A_{i}(x)A_{j}(0)\rangle-(-1)^{\tau/a_{\tau}}\langle B_{i}(x)B_{j}(0)\rangle,$
$\displaystyle A_{i}$ $\displaystyle=$
$\displaystyle\bar{\psi}\gamma_{i}\psi,\qquad\quad
B_{i}=\bar{\psi}\gamma_{4}\gamma_{5}\gamma_{i}\psi.$ (48)

To resolve this staggering, the spectral reconstruction is performed on even
and on odd time slices independently, obtaining two spectral functions,
$\rho^{\rm even,odd}_{ij}(\omega)$. In the sum

$\rho_{ij}(\omega)=\frac{1}{2}\left[\rho^{\rm even}_{ij}(\omega)+\rho^{\rm
odd}_{ij}(\omega)\right],$ (49)

the oscillating contribution cancels, and one can proceed with
$\rho_{ij}(\omega)$ as usual. This procedure limits, however, the number of
time slices effectively available by a factor of two. This issue does not
arise with Wilson-type fermions, which are hence the preferred choice.

  * •

local/conserved current: the simplest choice for the current operator is the
local one, $j^{\rm loc}_{x,i}=\bar{\psi}_{x}\gamma_{i}\psi_{x}$, with the
quark fields residing at the same lattice point $x$. This operator requires
renormalisation, i.e. the determination of a renormalisation factor $Z_{V}$.
The earliest contributions Gupta:2003zh ; Aarts:2007wj  did not determine this
renormalisation factor, which was accounted for in the systematic uncertainty.
Later studies using the local current did renormalise it properly.

The conserved current, of the point-split form (the expression below is for
Wilson fermions, $U_{x,i}$ is the gauge link in the spatial direction)

$j^{\rm
cons}_{x,i}=\frac{1}{2}\bar{\psi}_{x+\hat{\imath}}(1+\gamma_{i})U^{\dagger}_{x,i}\psi_{x}-\frac{1}{2}\bar{\psi}_{x}(1-\gamma_{i})U_{x,i}\psi_{x+\hat{\imath}},$
(50)

is the Noether current corresponding to a global phase symmetry of the fermion
lattice action and hence does not require renormalisation, even at finite
lattice spacing. For the current-current correlator, $\langle
j_{x,i}j_{y,i}\rangle$, combining two conserved currents, as in Refs.
Amato:2013naa ; Aarts:2014nba ; Astrakhantsev:2019zkr , eliminates the need to
compute the $Z_{V}$ factor, but it is also the most expensive numerically, due
to the need to invert more combinations of quark propagators. Ref.
Brandt:2015aqk  employs a mixed combination, $\langle j^{\rm cons}_{x,i}j^{\rm
loc}_{y,i}\rangle$, which is cheaper to evaluate and exactly conserved on one
side. It still requires knowledge of the $Z_{V}$ factor, which can e.g. be
obtained from $\langle j^{\rm loc}_{x,i}j^{\rm loc}_{y,i}\rangle/\langle
j^{\rm cons}_{x,k}j^{\rm loc}_{y,k}\rangle$.

  * •

the final column lists the methods, discussed above, employed to reconstruct
the spectral function and extract the conductivity. Ansätze and MEM have
traditionally been the most popular ones, with the Backus-Gilbert method and
Tikhonov regularisation being applied to this problem more recently.

Before moving to the results obtained so far, we note that none of the papers
listed above include the so-called disconnected contributions, discussed in
Sec. 2. The neglect of the disconnected contribution can be motivated in a
number of ways:

  * •

the contribution vanishes at very high $T$, since it is ${\cal
O}(\alpha_{s}^{3})$ at leading order in perturbation theory;

  * •

the contribution vanishes for three degenerate flavours, due to the sum over
the charges;

  * •

the contribution is expected to be noisy numerically.

This last comment is an excuse, rather than a motivation, and indeed, it would
be of interest to estimate the level of statistics required to compute a
signal in the $N_{f}=2$ or the non-degenerate $N_{f}=2+1$ case and verify e.g.
the first remark at high temperature. We also note that all studies discussed
here have been carried out at zero spatial momentum; the extension to nonzero
momentum is straightforward (see e.g. Ref. Aarts:2006cq ) and can provide an
additional handle on hydrodynamic behaviour at small $\omega$ and
$|{\mathbf{p}}|$.

After this overview of the lattice details, we are now in a position to
compare the conductivities computed in the references listed above. The
conductivity is normalised with the temperature (to make it dimensionless) and
with $C_{\rm em}=\sum_{f}(eq_{f})^{2}$, the sum over the charges squared, see
Eq. (16). The latter division allows one to compare e.g. the $N_{f}=2$ and
$2+1$ cases. We separately discuss the quenched results and the results with
dynamical quarks.

Figure 2: Results for the electrical conductivity, normalised as $\sigma/(TC_{\rm em})$, in quenched QCD ($N_{f}=0$) as a function of $T/T_{c}$. The early result Gupta:2003zh , $\sigma/(TC_{\rm em})\approx 7$, is not shown for clarity. Table 4: Estimates of the pseudocritical temperatures for the studies with dynamical quarks. ref. | $N_{f}$ |  $m_{\pi}$ [MeV] |  $T_{\rm pc}$ [MeV]
---|---|---|---
Brandt:2012jc ; Brandt:2015aqk  | $2$ | $270$ | $\sim 203$
Amato:2013naa ; Aarts:2014nba  | $2+1$ | $384(4)$ | $185(4)$
Astrakhantsev:2019zkr  | $2+1$ | $134.2(6)$ | $155(2)(3)$

Results in quenched QCD are shown in Fig. 2, as a function of $T/T_{c}$. The
very early result from Ref. Gupta:2003zh , $\sigma/(TC_{\rm em})\approx 7$, is
not included, since it is about a factor of $20$ larger than the other
results. The remaining four quenched studies are in good agreement, with a
value of $\sigma/(TC_{\rm em})\approx 0.2-0.5$. Note that Ref. Ding:2010ga
provides both a precise result, $\sigma/(TC_{\rm em})=0.37(1)$ at
$T/T_{c}=1.45$, and a more conservative range, indicated with the tallest
vertical green column. Although these are four studies, we note that they
emerge from two groups only, Ref. Aarts:2007wj  on the one hand and Refs.
Ding:2010ga ; Francis:2011bt ; Ding:2016hua  on the other hand, making the
agreement is perhaps less surprising. In any case, it is interesting that the
early quenched result of Ref. Aarts:2007wj , obtained using staggered quarks
without taking a continuum limit, remains to be consistent with the
renormalised continuum-extrapolated Wilson-clover results of Refs. Ding:2010ga
; Francis:2011bt ; Ding:2016hua . A second observation of interest is that
there appears to be very little temperature dependence in the temperature
range investigated, $1.1<T/T_{c}<3$. We remind the reader that in quenched QCD
the deconfinement transition is first-order, signalled by the spontaneous
breaking of the centre symmetry.

Figure 3: Results for electrical conductivity, normalised as $\sigma/(TC_{\rm
em})$, in QCD with $N_{f}=2$ and $2+1$ dynamical flavours as a function of
temperature in MeV.

We now turn to the dynamical results, shown in Fig. 3 as a function of
temperature in MeV. In this case, the thermal transition is a smooth
crossover. It should be noted that the crossover temperature is slightly
different between the various studies, due to the difference in the pion mass
and the number of flavours, see Table 4 (note that Ref. Astrakhantsev:2019zkr
follows the lattice formulation and choice of parameters of Refs. Aoki:2009sc
; Borsanyi:2010cj . Moreover, Ref. Astrakhantsev:2019zkr  has results at two
lattice spacings; the red star symbols denote the results at the smaller
spacing, $N_{\tau}=16$). In particular, the transition temperatures in Refs.
Amato:2013naa ; Aarts:2014nba ; Brandt:2015aqk  are higher than in Ref.
Astrakhantsev:2019zkr , with simulations in the latter being at the physical
point. In Fig. 3 a temperature-dependent $\sigma/T$ can be observed, with a
reduction in the vicinity of the thermal crossover, in contrast to the
quenched case. This temperature dependence is especially visible in the data
from Refs. Amato:2013naa ; Aarts:2014nba  and Ref. Brandt:2015aqk , which have
results at a number of temperatures (8 and 4 values respectively). The
difference in crossover temperatures might explain the slightly lower values
for the conductivity at $T=200$ MeV in Refs. Amato:2013naa ; Aarts:2014nba ;
Brandt:2015aqk , compared to Ref. Astrakhantsev:2019zkr , as in the former
cases 200 MeV is closer to the pseudocritical temperatures and the thermal
crossover region. As mentioned above, in quenched QCD, with a first-order
transition, no reduction of the conductivity above the critical temperature is
seen. These observations suggest that the reduction of the conductivity is due
to the smooth transition to the hadronic phase. We also note that Refs.
Brandt:2012jc ; Brandt:2015aqk  are from the same group; hence the data point
indicated by the square is effectively superseded by those indicated with the
circles. This allows us to conclude that there is good consistency between the
various studies, which is a nontrivial result, given the difference in lattice
formulations, lattice geometries and inversion methods employed. Finally we
observe that at the highest temperatures studied the value for
$\sigma/(TC_{\rm em})$ is comparable to the one obtained in the quenched case.

In order to judge whether the observed magnitude of the conductivity signifies
strong- or weak-coupling behaviour, we note that at weak coupling (including
only QCD processes) the usual expectation is that $\sigma/(TC_{\rm em})\sim
1/g^{4}\ln 1/g$ Arnold:2000dr , which is much larger than 1 in the region
where the weak-coupling analysis is valid, namely at asymptotically high
temperatures. A useful benchmark at strong coupling comes from holography, for
the charge diffusion coefficient $D=\sigma/\chi_{Q}$, where $\chi_{Q}$ is the
charge susceptibility. The characteristic result at strong coupling is
$D=1/(2\pi T)$ in ${\cal N}=4$ Yang- Mills theory at nonzero temperature
Policastro:2002se ; Teaney:2006nc . In Ref. Aarts:2014nba  the temperature
dependence of $D$ was computed in a self-contained manner, i.e. by also
computing $\chi_{Q}$ within the same lattice QCD setup, with the result that
$0.5<2\pi TD<2$, compatible with the holographic order of magnitude at strong
coupling. Moreover, it was observed that $2\pi TD$ has a minimum in the
crossover region, see Fig. 14 of Ref. Aarts:2014nba .

Before concluding this section, we note that the lattice data of Ref.
Ding:2010ga  have been re-analysed in two papers. Ref. Burnier:2012ts
employed the approach of Ref. Cuniberti , see also Ref. Burnier:2011jq , in
which short-distance divergences are subtracted from the Euclidean correlator.
Additional insight on the ultraviolet asymptotics of the thermal contribution
to the spectral function is taken from Ref. CaronHuot:2009ns , such that only
the contribution of the vacuum spectral function needs to be subtracted. For
this, a 5-loop computation of the vector current correlator in vacuum is
employed. A smaller result for the conductivity and diffusion coefficient are
found with respect to Ref. Ding:2010ga , namely, at $T/T_{c}=1.45$,

$\sigma/(TC_{\rm em})\gtrsim~{}0.1,\qquad\quad 2\pi TD\gtrsim 0.8,$ (51)

which is indeed smaller by a factor of 3 for the conductivity. It is stated
Burnier:2012ts  that the results in Eq. (51) should be interpreted as lower
bounds.

The data of Refs. Ding:2010ga  and Brandt:2015aqk  has also been re-analysed
in Refs. Gubler:2016hnf  and Gubler:2017qbs  respectively, using thermal sum
rules to constrain the Ansätze used in the fits. In Ref. Gubler:2016hnf , a
higher result was found for quenched study at $T/T_{c}=1.45$, namely

$\sigma/(TC_{\rm em})\sim 0.57.$ (52)

In Ref. Gubler:2017qbs  the $N_{f}=2$ data given in Ref. Brandt:2015aqk  at
four temperatures was re-analysed. Approximate agreement was found, with the
important caveat that the fits were seen not to be as stable as desired.

Even though good consistency between the various studies can be seen,
nevertheless continuing uncertainty in spectral reconstruction remains, with
an ongoing need to further develop methods for analytical continuation,
emphasising robustness and quantification of underlying uncertainties.

##  5 External conditions

The electrical conductivity, as well as other transport coefficients, may be
studied under other external conditions than temperature, such as in the
presence of an external magnetic field or at nonzero quark density. In this
section we briefly mention some related developments in nonabelian gauge
theories.

The first attempts to investigate the dependence of the electrical
conductivity on an external magnetic field using lattice simulations were
performed in the quenched SU(2) theory in Refs. Buividovich:2010tn ;
Buividovich:2010qe , using the overlap Dirac operator with exact chiral
symmetry in the current-current correlator. MEM is used for spectral
reconstruction. The emphasis is on the conductivity in the presence of an
external magnetic field, both in the confined and the deconfined phase, and on
the quark mass dependence. It is found that in the confined phase the external
magnetic field induces a nonzero electric conductivity along the direction of
the field, while in the deconfined phase no sizable dependence on the magnetic
field is observed.

In Ref. Buividovich:2020el  the conductivity is studied in the SU(2) gauge
theory with dynamical quarks at nonzero density, which is feasible due the
absence of a sign problem in this theory. Gauge configurations are generated
with dynamical staggered quarks, while current-current correlators are
computed with Wilson-Dirac and Domain Wall fermions, tuned in such a way to
match the pion mass of the ensembles. The conductivity is extracted via
several methods, including the Backus-Gilbert method with the use of Tikhonov
regularisation. At small quark chemical potential $\mu$, the dependence on
chemical potential is considered via the expansion

$\frac{\sigma(\mu)}{\sigma(0)}=1+c(T)\frac{\mu^{2}}{T^{2}}+{\cal
O}\left(\frac{\mu^{4}}{T^{4}}\right),$ (53)

and it is found that the maximal value of the second-order coefficient,
$c(T)\approx 0.15(5)$, is reached in the vicinity of the chiral crossover.
Hence the coefficient $c(T)$ is quite small, and even at $\mu/T\approx 1$ the
conductivity changes no more than 15–20% compared to its zero-density value.
As for the large density region, QCD with $N_{c}=2$ and $3$ colours differ,
due to the formation of a diquark condensate in the former. Nevertheless, at
smaller $\mu$ the SU(2) theory can provide qualitative insights for real QCD.

Also in QCD (with $N_{c}=3$ and $N_{f}=2+1$) the electrical conductivity has
been studied in the presence of an external constant uniform magnetic field
Astrakhantsev:2019zkr . The ${\mathbf{B}}={\mathbf{0}}$ results of this
reference have been discussed above; with nonzero ${\mathbf{B}}$ field, it is
found that the conductivity rises in the direction parallel to the magnetic
field and decreases in the transverse direction. This may potentially be
explained by the Chiral Magnetic Effect Kharzeev:2015znc  and
magnetoresistance.

##  6 Conclusion

In this paper we reviewed the status of the electrical conductivity in the
quark-gluon plasma, as seen through nonperturbative lattice QCD simulations.
After an overview of basic definitions and expectations, we listed several
methods that have been used for spectral reconstruction, the main challenge in
this endeavour. No method has yet reached full acceptance, due to the apparent
lack of robustness and handle on systematic uncertainties. This remains
therefore the outstanding challenge to be tackled. It is also noted that none
of the results include the disconnected contributions yet, for reasons
discussed in Sec. 4. It would be worthwhile to estimate the importance of
those eventually.

Nevertheless, a comparison between the existing lattice studies, presented in
Sec. 4, reveals a noticeable consistency, which is encouraging, given the
difference in lattice formulations, lattice geometries (in particular the
number of temporal points), and reconstruction methods employed. Taking the
results at face value, the main findings are

  * •

in quenched ($N_{f}=0$) QCD $\sigma/T$ appears to have very little temperature
dependence in the temperature range investigated, $1.1<T/T_{c}<3$. The
magnitude is approximately $0.2\lesssim\sigma/(TC_{\rm em})\lesssim 0.5$,
where $C_{\rm em}$ is the sum over electric charges squared appearing in the
electromagnetic current-current correlator;

  * •

in QCD with $N_{f}=2$ and $2+1$ dynamical flavours, the main finding is a
noticeable reduction of $\sigma/T$ in the vicinity of the thermal crossover,
compared to its value at higher temperatures in the QGP. This should be
contrasted with the quenched case. This effect has been observed by two groups
independently and is further (indirectly) supported by simulations at the
physical point by a third group. One possible interpretation is that the
reduction of the conductivity is due to the smooth transition to the hadronic
phase. This might be of interest for phenomenology. It is further noted that
at the highest temperatures studied the value for $\sigma/(TC_{\rm em})$ is
comparable to the one obtained in the quenched case. Overall, the magnitude of
the conductivity is compatible with the plasma being strongly coupled, using
the comparison of the charge diffusion coefficient $D=\sigma/\chi_{Q}$, where
$\chi_{Q}$ is the charge susceptibility, with the one obtained in holography.

So far most studies have focused on the quark-gluon plasma and the crossover
region. Deeper in the hadronic phase the conductivity should be dominated by
the lightest charged hadrons. So far, lattice studies have not given a
detailed study of this regime, possibly because the signal is hard to detect.
We refer to Ref. FernandezFraile:2009mi  for an overview from the perspective
of chiral perturbation theory. Studies in the presence of an external magnetic
field or at finite density require further attention. While direct access to
the latter is not feasible in QCD due to the sign problem, a Taylor series
expansion in the powers of $\mu/T$ is possible, although it is expected to be
noisy numerically Aarts:2020lvh . Another interesting possibility is the
analysis of the current-current correlator at imaginary chemical potential.

###### Acknowledgements.

We thank Mikko Laine and Pavel Buividovich for discussion. GA is grateful for
collaboration with Chris Allton, Ale Amato, Justin Foley, Pietro Giudice,
Simon Hands, Seyong Kim, and Jon-Ivar Skullerud on some of the results
reviewed here. This work is supported by STFC via grant ST/P00055X/1, the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme via grant agreement No. 813942, and by COST
Action CA15213 THOR. A.N. also acknowledges the support by RFBR grant No.
18-02-40126.

##  Appendix A Green’s functions

For the convenience of the reader we collect in this Appendix some relations
between the various two-point functions for an operator
$O(x)=O(t_{x},{\mathbf{x}})$. These relations are well known LeBellac ;
Laine:2016hma .

Let us start with the retarded and advanced Green’s functions

$G_{R}(x-y)=i\theta(t_{x}-t_{y})\langle[O(x),O^{\dagger}(y)]\rangle=G_{A}(y-x),$
(54)

and the spectral function

$\displaystyle\rho(x-y)$ $\displaystyle=$
$\displaystyle\langle[O(x),O^{\dagger}(y)]\rangle$ (55) $\displaystyle=$
$\displaystyle-i\left[G_{R}(x-y)-G_{A}(x-y)\right].$

Expectation values are taken in thermal equilibrium, which explains the $x-y$
dependence. After going to momentum space and using the identity

$\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\frac{e^{-i\omega
t}}{\omega+i\epsilon}=-i\theta(t),$ (56)

we arrive at the dispersion relation

$G_{R}(\omega,{\mathbf{p}})=\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\,\frac{\rho(\omega^{\prime},{\mathbf{p}})}{\omega^{\prime}-\omega-i\epsilon}.$
(57)

Employing the identity

$\frac{1}{x+i\epsilon}-\frac{1}{x-i\epsilon}=\frac{-2i\epsilon}{x^{2}+\epsilon^{2}}\to-2i\pi\delta(x),$
(58)

then yields the important relation,

$\rho(\omega,{\mathbf{p}})=-i\left[G_{R}(\omega,{\mathbf{p}})-G_{A}(\omega,{\mathbf{p}})\right]=2{\mathrm{Im}}\,G_{R}(\omega,{\mathbf{p}}),$
(59)

i.e. the spectral function is twice the imaginary part of retarded Green
function, or equivalently the discontinuity across the real axis.

The Euclidean correlator,

$G_{E}(\tau,{\mathbf{x}})=\langle
O(\tau,{\mathbf{x}})O^{\dagger}(0,{\mathbf{0}})\rangle,$ (60)

with $0\leq\tau<1/T$, is written in momentum space as

$\displaystyle G_{E}(\omega_{n},{\mathbf{p}})$ $\displaystyle=$
$\displaystyle\int_{0}^{1/T}d\tau\,e^{i\omega_{n}\tau}G_{E}(\tau,{\mathbf{p}}),$
(61) $\displaystyle G_{E}(\tau,{\mathbf{p}})$ $\displaystyle=$ $\displaystyle
T\sum_{n}e^{-i\omega_{n}\tau}G_{E}(\omega_{n},{\mathbf{p}}),$ (62)

where $\omega_{n}=2\pi nT$, $n\in\mathbb{Z}$, are the Matsubara frequencies
(we consider bosonic operators here). By analyticity, it satisfies a similar
dispersion relation as above,

$G_{E}(\omega_{n},{\mathbf{p}})=\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\,\frac{\rho(\omega^{\prime},{\mathbf{p}})}{\omega^{\prime}-i\omega_{n}},$
(63)

leading to the important relation

$G_{R}(\omega,{\mathbf{p}})=G_{E}(i\omega_{n}\to w+i\epsilon,{\mathbf{p}}).$
(64)

If a Euclidean correlator is known analytically, the spectral function can be
obtained following the sequence

$G_{E}(\tau,{\mathbf{x}})\to G_{E}(\omega_{n},{\mathbf{p}})\to
G_{R}(\omega,{\mathbf{p}})\to\rho(\omega,{\mathbf{p}}).$ (65)

Unfortunately, this path is not accessible with numerically determined
correlators on a finite number of points in the temporal direction.

Instead we will relate the correlator and the spectral function via a Laplace
transform, generalised to nonzero temperature. Going back to Euclidean time,
we find, using Eq. (63),

$G(\tau,{\mathbf{p}})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,\tilde{K}(\tau,\omega)\rho(\omega,{\mathbf{p}}),$
(66)

with the kernel

$\tilde{K}(\tau,\omega)=T\sum_{n}\frac{e^{-i\omega_{n}\tau}}{\omega-i\omega_{n}}=e^{-\omega\tau}[1+n_{B}(\omega)].$
(67)

Here $n_{B}(\omega)=1/[\exp(\omega/T)-1]$ is the Bose-Einstein distribution
and we have taken $0<\tau<1/T$.

For hermitian operators $O^{\dagger}(x)=O(x)$, the spectral function is odd in
$\omega$,

$\rho(-\omega,{\mathbf{p}})=-\rho(\omega,{\mathbf{p}}).$ (68)

Hence in Eq. (66) only the odd part of the kernel $\tilde{K}(\tau,\omega)$
survives, and we arrive at the standard integral relation

$G(\tau,{\mathbf{p}})=\int_{0}^{\infty}\frac{d\omega}{2\pi}\,K(\tau,\omega)\rho(\omega,{\mathbf{p}}),$
(69)

with

$\displaystyle K(\tau,\omega)$ $\displaystyle=$
$\displaystyle\tilde{K}(\tau,\omega)-\tilde{K}(\tau,-\omega)$ (70)
$\displaystyle=$ $\displaystyle
e^{-\omega\tau}[1+n_{B}(\omega)]+e^{\omega\tau}n_{B}(\omega)$ $\displaystyle=$
$\displaystyle\frac{\cosh[\omega(\tau-1/2T)]}{\sinh(\omega/2T)},$ (71)

where we used the identity

$n_{B}(\omega)+n_{B}(-\omega)+1=0.$ (72)

## References

  * (1) W. T. Deng and X. G. Huang, Phys. Rev. C 85 (2012) 044907 [arXiv:1201.5108 [nucl-th]].
  * (2) Y. Hirono, M. Hongo and T. Hirano, Phys. Rev. C 90 (2014) no.2, 021903 [arXiv:1211.1114 [nucl-th]].
  * (3) K. Tuchin, Adv. High Energy Phys. 2013 (2013) 490495 [arXiv:1301.0099 [hep-ph]].
  * (4) Y. Yin, Phys. Rev. C 90 (2014) no.4, 044903 [arXiv:1312.4434 [nucl-th]].
  * (5) D. E. Kharzeev, J. Liao, S. A. Voloshin and G. Wang, Prog. Part. Nucl. Phys. 88 (2016) 1 [arXiv:1511.04050 [hep-ph]].
  * (6) S. Pratt and C. Plumberg, arXiv:1904.11459 [nucl-th].
  * (7) G. Aarts and J. M. Martínez Resco, JHEP 0211 (2002) 022 [hep-ph/0209048].
  * (8) J. S. Gagnon and S. Jeon, Phys. Rev. D 75 (2007) 025014 Erratum: [Phys. Rev. D 76 (2007) 089902] [hep-ph/0610235].
  * (9) P. B. Arnold, G. D. Moore and L. G. Yaffe, JHEP 0011 (2000) 001 [hep-ph/0010177].
  * (10) P. B. Arnold, G. D. Moore and L. G. Yaffe, JHEP 0305 (2003) 051 [hep-ph/0302165].
  * (11) W. Cassing, O. Linnyk, T. Steinert and V. Ozvenchuk, Phys. Rev. Lett. 110 (2013) no.18, 182301 [arXiv:1302.0906 [hep-ph]].
  * (12) R. Marty, E. Bratkovskaya, W. Cassing, J. Aichelin and H. Berrehrah, Phys. Rev. C 88 (2013) 045204 [arXiv:1305.7180 [hep-ph]].
  * (13) G. Policastro, D. T. Son and A. O. Starinets, JHEP 0209 (2002) 043 [hep-th/0205052].
  * (14) D. Teaney, Phys. Rev. D 74 (2006) 045025 [hep-ph/0602044].
  * (15) S. Gupta, Phys. Lett. B 597 (2004) 57 [hep-lat/0301006].
  * (16) G. Aarts, C. Allton, J. Foley, S. Hands and S. Kim, Phys. Rev. Lett. 99 (2007) 022002 [hep-lat/0703008 [HEP-LAT]].
  * (17) H.-T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann and W. Soeldner, Phys. Rev. D 83 (2011) 034504 [arXiv:1012.4963 [hep-lat]].
  * (18) A. Francis and O. Kaczmarek, Prog. Part. Nucl. Phys. 67 (2012) 212 [arXiv:1112.4802 [hep-lat]].
  * (19) B. B. Brandt, A. Francis, H. B. Meyer and H. Wittig, JHEP 1303 (2013) 100 [arXiv:1212.4200 [hep-lat]].
  * (20) A. Amato, G. Aarts, C. Allton, P. Giudice, S. Hands and J. I. Skullerud, Phys. Rev. Lett. 111 (2013) no.17, 172001 [arXiv:1307.6763 [hep-lat]].
  * (21) G. Aarts, C. Allton, A. Amato, P. Giudice, S. Hands and J. I. Skullerud, JHEP 1502 (2015) 186 [arXiv:1412.6411 [hep-lat]].
  * (22) B. B. Brandt, A. Francis, B. Jäger and H. B. Meyer, Phys. Rev. D 93 (2016) no.5, 054510 [arXiv:1512.07249 [hep-lat]].
  * (23) H. T. Ding, O. Kaczmarek and F. Meyer, Phys. Rev. D 94 (2016) no.3, 034504 [arXiv:1604.06712 [hep-lat]].
  * (24) N. Y. Astrakhantsev, V. V. Braguta, M. D’Elia, A. Y. Kotov, A. A. Nikolaev and F. Sanfilippo, arXiv:1910.08516 [hep-lat].
  * (25) H. B. Meyer, New J. Phys. 13 (2011) 035008 [arXiv:1012.0234 [hep-lat]].
  * (26) D. Banerjee, S. Datta, R. Gavai and P. Majumdar, Phys. Rev. D 85 (2012) 014510 [arXiv:1109.5738 [hep-lat]].
  * (27) A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus and H. Ohno, Phys. Rev. D 92 (2015) no.11, 116003 [arXiv:1508.04543 [hep-lat]].
  * (28) H. B. Meyer, Phys. Rev. D 76 (2007) 101701 [arXiv:0704.1801 [hep-lat]].
  * (29) N. Astrakhantsev, V. Braguta and A. Kotov, JHEP 1704 (2017) 101 [arXiv:1701.02266 [hep-lat]].
  * (30) S. Borsányi et al., Phys. Rev. D 98 (2018) no.1, 014512 [arXiv:1802.07718 [hep-lat]].
  * (31) H. B. Meyer, Phys. Rev. Lett. 100 (2008) 162001 [arXiv:0710.3717 [hep-lat]].
  * (32) N. Y. Astrakhantsev, V. V. Braguta and A. Y. Kotov, Phys. Rev. D 98 (2018) no.5, 054515 [arXiv:1804.02382 [hep-lat]].
  * (33) G. Aarts and J. M. Martínez Resco, Nucl. Phys. B 726 (2005) 93 [hep-lat/0507004].
  * (34) G. Aarts and J. M. Martínez Resco, JHEP 0204 (2002) 053 [hep-ph/0203177].
  * (35) G. Aarts and J. M. Martínez Resco, Nucl. Phys. Proc. Suppl. 119 (2003) 505 [hep-lat/0209033].
  * (36) P. Petreczky and D. Teaney, Phys. Rev. D 73 (2006), 014508 [arXiv:hep-ph/0507318 [hep-ph]].
  * (37) H. B. Meyer, Eur. Phys. J. A 47 (2011), 86 [arXiv:1104.3708 [hep-lat]].
  * (38) S. Datta, F. Karsch, P. Petreczky and I. Wetzorke, Phys. Rev. D 69 (2004), 094507 [arXiv:hep-lat/0312037 [hep-lat]].
  * (39) G. Aarts, C. Allton, M. B. Oktay, M. Peardon and J. I. Skullerud, Phys. Rev. D 76 (2007), 094513 [arXiv:0705.2198 [hep-lat]].
  * (40) S. Caron-Huot, Phys. Rev. D 79 (2009) 125009 [arXiv:0903.3958 [hep-ph]].
  * (41) P. Gubler and D. Satow, Phys. Rev. D 94 (2016) no.9, 094042 [arXiv:1602.08265 [hep-ph]].
  * (42) P. Gubler and D. Satow, Phys. Rev. D 96 (2017) no.11, 114028 [arXiv:1710.02256 [hep-ph]].
  * (43) F. Karsch, E. Laermann, P. Petreczky and S. Stickan, Phys. Rev. D 68 (2003), 014504 [arXiv:hep-lat/0303017 [hep-lat]].
  * (44) Y. Burnier and M. Laine, Eur. Phys. J. C 72 (2012) 1902 [arXiv:1201.1994 [hep-lat]].
  * (45) R. K. Bryan, Eur. Biophys. J. 18 (1990) 165.
  * (46) M. Asakawa, T. Hatsuda and Y. Nakahara, Prog. Part. Nucl. Phys. 46 (2001), 459-508 [arXiv:hep-lat/0011040 [hep-lat]].
  * (47) A. Tikhonov, Soviet Math. Dokl. 5 (1963) 1035/1038 [http://ci.nii.ac.jp/naid/10004315593/en/].
  * (48) G. Cuniberti, E. De Micheli and G.A. Viano, Commun. Math. Phys. 216 (2001) 59 [arXiv:cond-mat/0109175].
  * (49) Y. Burnier, M. Laine and L. Mether, Eur. Phys. J. C 71 (2011) 1619 [arXiv:1101.5534 [hep-lat]].
  * (50) Y. Burnier and A. Rothkopf, Phys. Rev. Lett. 111 (2013), 182003 [arXiv:1307.6106 [hep-lat]].
  * (51) A. Rothkopf, J. Comput. Phys. 238 (2013), 106-114 [arXiv:1110.6285 [physics.comp-ph]].
  * (52) R. A. Tripolt, P. Gubler, M. Ulybyshev and L. Von Smekal, Comput. Phys. Commun. 237 (2019) 129 [arXiv:1801.10348 [hep-ph]].
  * (53) L. Kades, J. M. Pawlowski, A. Rothkopf, M. Scherzer, J. M. Urban, S. J. Wetzel, N. Wink and F. Ziegler, (2019) [arXiv:1905.04305].
  * (54) L.-F. Arsenault, R. Neuberg, L.A. Hannah and A.J. Millis, Inverse Problems 33 (2017) 115007 [arXiv:1612.04895 [cond-mat]].
  * (55) S. Offler, G. Aarts, C. Allton, J. Glesaaen, B. Jäger, S. Kim, M. P. Lombardo, S. M. Ryan and J. I. Skullerud, PoS LATTICE2019 (2019) 076 [arXiv:1912.12900 [hep-lat]].
  * (56) G. Aarts, C. Allton, J. Foley, S. Hands and S. Kim, PoS LAT2006 (2006), 134 [arXiv:hep-lat/0610061 [hep-lat]].
  * (57) Y. Aoki, S. Borsányi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 06 (2009), 088 [arXiv:0903.4155 [hep-lat]].
  * (58) S. Borsányi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, JHEP 1011 (2010) 077 [arXiv:1007.2580 [hep-lat]].
  * (59) P. V. Buividovich, M. N. Chernodub, D. E. Kharzeev, T. Kalaydzhyan, E. V. Luschevskaya and M. I. Polikarpov, Phys. Rev. Lett. 105 (2010) 132001 [arXiv:1003.2180 [hep-lat]].
  * (60) P. V. Buividovich and M. I. Polikarpov, Phys. Rev. D 83 (2011) 094508 [arXiv:1011.3001 [hep-lat]].
  * (61) P. V. Buividovich, D. Smith and L. von Smekal, arXiv:2007.05639 [hep-lat].
  * (62) D. Fernandez-Fraile and A. Gomez Nicola, Eur. Phys. J. C 62 (2009), 37-54 [arXiv:0902.4829 [hep-ph]].
  * (63) A. A. Nikolaev, G. Aarts, C. Allton, D. De Boni, J. Glesaaen, S. Hands, B. Jäger, J. I. Skullerud and L. K. Wu, PoS LATTICE2019 (2019) 077 [arXiv:2001.04415 [hep-lat]].
  * (64) M. Le Bellac, “Thermal Field Theory,” Cambridge UNiversity Press (1996).
  * (65) M. Laine and A. Vuorinen, Lect. Notes Phys. 925 (2016) pp.1 [arXiv:1701.01554 [hep-ph]].