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.  [1,2,3,4,5,6] 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 [7, 8] and kinetic theory in QCD [9, 10] or effective models [11, 12] to holographic methods at strong coupling [13, 14]. 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 \({{{\mathcal {O}}}}(10)\) papers which have attempted to compute the conductivity on the lattice [15,16,17,18,19,20,21,22,23,24]. 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 Sect. 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 [17,18,19, 22, 23], sum rule constraints [19, 25, 26], Bayesian approaches such as the Maximum Entropy Method [15,16,17,18, 20, 21], the Backus-Gilbert method [22, 24], and analytical continuation of data after short-distance vacuum subtractions [27,28,29]. 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 [30,31,32] and neither other transport coefficients such as the shear [33,34,35] and bulk viscosities [36, 37].

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 Sect.  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 Sect.  4, including a comparison between the values of \(\sigma \) obtained so far. Some related developments are summarised in Sect. 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

$$\begin{aligned} j_\mu ^{\mathrm{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).\nonumber \\ \end{aligned}$$
(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 ^{\mathrm{em}}}^\dagger (x)=j_\mu ^{\mathrm{em }}(x)\).

The electrical conductivity \(\sigma \) indicates the linear relationship between the current density and an electric field, \(j_i^{\mathrm{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. [9]. More precisely, the conductivity is proportional to the slope of the current-current spectral function in thermal equilibrium,

$$\begin{aligned} \rho _{\mu \nu }^{\mathrm{em}}(\omega ,{{\mathbf {p}}}) = \int d^4x\, e^{i\omega t-{{\mathbf {p}}}\cdot {{\mathbf {x}}}} \langle [j_\mu ^{\mathrm{em}}(t,{{\mathbf {x}}}), j_\nu ^{\mathrm{em}}(0,{{\mathbf {0}}})]\rangle , \end{aligned}$$
(2)

at vanishing energy and momentum, i.e.,

$$\begin{aligned} \sigma = \frac{1}{6}\frac{\partial }{\partial \omega }\rho ^{\mathrm{em}}_{ii}(\omega ,{{\mathbf {0}}})\Big |_{\omega =0}. \end{aligned}$$
(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 [9] to the charge diffusion coefficient D, according to the Einstein relation,

$$\begin{aligned} \sigma = \chi _Q D, \end{aligned}$$
(4)

where \(\chi _Q\) is the charge susceptibility,

$$\begin{aligned} \chi _Q = \frac{1}{TV} \left\langle (Q-\langle Q\rangle )^2\right\rangle . \end{aligned}$$
(5)

Here V is the spatial volume and Q is the total charge, i.e. the volume integral of \(j_0^{\mathrm{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,

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

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

$$\begin{aligned} G_{\mu \nu }^{\mathrm{em}}(\tau ,{{\mathbf {p}}}) = \int _0^\infty \frac{d\omega }{2\pi }\, K(\tau ,\omega )\rho _{\mu \nu }^{\mathrm{em}}(\omega , {{\mathbf {p}}}), \end{aligned}$$
(7)

with the kernel

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

Here the Euclidean time \(0\le \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.

Lattice QCD simulations include quarks scattering with gluons, but the electromagnetic field and other charge carriers (leptons) are not included. The QCD computation will hence yield a result valid to all orders in \(\alpha _s\) and to leading order in \(\alpha _{\mathrm{em}}\), in principle (see, however, the discussion on disconnected contributions further down in this Section). Since the conductivity reviewed here includes the contribution from the strong interaction, it yields insight into the strongly coupled nature of the quark-gluon plasma.

Fig. 1
figure 1

One-loop contribution to the current-current spectral function

Full size image

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 an expansion in \(\alpha _s\) (see Fig. 1), the spectral function, for a single flavour with mass m, reads [38]Footnote 1

$$\begin{aligned} \rho _{ii}(\omega )= & {} 2\pi N_c I\omega \delta (\omega ) \nonumber \\&+ \frac{N_c}{2\pi } \theta (\omega ^2-4m^2)\sqrt{\frac{\omega ^2-4m^2}{\omega ^2}} \nonumber \\&\times \left( \omega ^2+2m^2\right) \left[ 1-2n_F\left( \frac{\omega }{2}\right) \right] , \end{aligned}$$
(9)

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

$$\begin{aligned} I = -4 \int \frac{d^3k}{(2\pi )^3}\, n_F'(\omega _{{\mathbf {k}}}) \frac{k^2}{\omega _{{\mathbf {k}}}^2}, \end{aligned}$$
(10)

with \(\omega _{{\mathbf {k}}}=\sqrt{{{\mathbf {k}}}^2+m^2}\). For massless quarks this evaluates as

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

Since

$$\begin{aligned} 1-2n_F\left( \frac{\omega }{2}\right) = \tanh \left( \frac{\omega }{4T}\right) , \end{aligned}$$
(12)

the spectral function is odd, \(\rho _{ii}(-\omega )=-\rho _{ii}(\omega )\), as it should be. Below we take \(\omega \ge 0\). The corresponding Euclidean correlator reads, in the massless limit [38],

$$\begin{aligned} G_{ii}(\tau )= N_c T^3 \left[ \frac{1}{3}+ \frac{3u+u\cos (2u)-2\sin (2u)}{\sin ^3(u)} \right] , \end{aligned}$$
(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_cT^3(1/3+2/3)=N_cT^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 \rightarrow 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}}}\rightarrow {{\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,

$$\begin{aligned} 2\pi N_c I \omega \delta (\omega ) \rightarrow \rho _{\mathrm{trans}}(\omega ) = A_{\mathrm{trans}} \frac{\gamma \omega }{\omega ^2+\gamma ^2}, \end{aligned}$$
(14)

where \(\gamma \) is proportional to the collisional scattering width or the inverse mean free path, \(A_{\mathrm{trans}}(=2N_c I)\) is the overall coefficient, and \(\sigma \propto A_{\mathrm{trans}}/\gamma \). The origin of the transport peak can be seen in various ways, using e.g. Feynman diagrams and pinching poles [39, 40] or kinetic theory [41, 42]. 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 [39]. In the limit that \(\omega < \varLambda \ll T\), the kernel (8) simplifies to \(2T/\omega \), and integrating the transport peak (14) in Eq. (7) yields

$$\begin{aligned} G_{\mathrm{trans}}(\tau )= & {} \int _0^{\varLambda } \frac{d\omega }{2\pi }\, K(\tau ,\omega )\rho _{\mathrm{trans}}(\omega ) \nonumber \\&\sim A_{\mathrm{trans}} \int _0^{\varLambda } \frac{d\omega }{2\pi }\, \frac{2T}{\omega } \frac{\gamma \omega }{\omega ^2+\gamma ^2} \sim \frac{1}{2}A_{\mathrm{trans}} T,\nonumber \\ \end{aligned}$$
(15)

where in the last expression the cutoff \(\varLambda \) 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_{\mathrm{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 \({{{\mathcal {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

$$\begin{aligned} C_{\mathrm{em}} = \sum _f (e q_f)^2 = {\left\{ \begin{array}{ll} \frac{5}{9} e^2 \qquad (u, d), \\ \frac{6}{9} e^2 \qquad (u, d, s). \end{array}\right. } \end{aligned}$$
(16)

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

$$\begin{aligned} C_{\mathrm{em}}^{\mathrm{disc}} = \left( \sum _f e q_f \right) ^2 = {\left\{ \begin{array}{ll} \frac{1}{9} e^2 \qquad (u, d), \\ 0 \qquad \quad (u, d, s). \end{array}\right. } \end{aligned}$$
(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. [19]. 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.

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 _{\mathrm{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 \rightarrow 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 Sect. 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\le \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 [43, 44]. 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

$$\begin{aligned} G_{\mathrm{recon}}(\tau , T;T_0) = \int _0^\infty \frac{d\omega }{2\pi }\, K(\tau ,\omega ;T) \rho (\omega ;T_0). \end{aligned}$$
(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 [19] to constrain the current-current spectral function at nonzero temperature. Defining the difference between the finite and zero-temperature spectral function as

$$\begin{aligned} \varDelta \rho (\omega ,T) \equiv \rho _{ii}(\omega ,T) - \rho _{ii}(\omega ,0), \end{aligned}$$
(19)

one finds the sum rule [19], in the thermodynamic limit,

$$\begin{aligned} \int _{-\infty }^\infty \frac{d\omega }{\omega } \varDelta \rho (\omega ,T)=0. \end{aligned}$$
(20)

Note that the zero-temperature \(\omega ^2\) contribution cancels in the subtraction (19). Since the Operator Product Expansion (OPE) predicts [45] 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. [25, 26].

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 [38, 46];

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

References [17, 18, 23] employ an Ansatz combining a transport peak and the expected perturbative continuum behaviour (for massless quarks) in the deconfined phase,

$$\begin{aligned} \rho (\omega ) = \rho _{\mathrm{trans}}(\omega ) + \rho _{\mathrm{pert}}(\omega ), \end{aligned}$$
(21)

where \(\rho _{\mathrm{trans}}(\omega )\) is given in Eq. (14) and

$$\begin{aligned} \rho _{\mathrm{pert}}(\omega ) = \frac{3}{2\pi } A_{\mathrm{pert}} \omega ^2 \left[ 1-2n_F\left( \frac{\omega }{2}\right) \right] , \end{aligned}$$
(22)

c.f. Eq. (9). The three temperature-dependent parameters are the coefficients \(A_{\mathrm{trans, pert}}\) and the width \(\gamma \) of the transport peak. Note that \(A_{\mathrm{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 [14], Ref. [23] 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 \) [23].

Reference [19] employs a related Ansatz for the subtracted spectral function, \(\varDelta \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

$$\begin{aligned} \varDelta \rho (\omega ) = \rho _{\mathrm{trans}}(\omega ) + \varDelta \rho _{\mathrm{pert}}(\omega ) + \rho _{\mathrm{bound}}(\omega ), \end{aligned}$$
(23)

with the new term

$$\begin{aligned} \rho _{\mathrm{bound}}(\omega ) = A_{\mathrm{bound}}\frac{2 g_B \tanh (\omega /T)^3}{4(\omega -m_B)^2 + g_B^2}. \end{aligned}$$
(24)

Here \(m_B, g_B\) and \(A_{\mathrm{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 \rightarrow 0\) limit and decays as \(1/\omega ^2\) at large \(\omega \), as predicted by the OPE [45]. It is also noted in Refs.  [19, 29] that the transport peak (14) in fact violates this condition. Hence it is proposed [19] to modify it as

$$\begin{aligned} \rho _{\mathrm{trans, mod}}(\omega ) = A_{\mathrm{trans}} \frac{T\tanh (\omega /T) \gamma }{\omega ^2+\gamma ^2}\,, \end{aligned}$$
(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\) term,

$$\begin{aligned} \varDelta \rho _{\mathrm{pert}}(\omega ) = \rho _{\mathrm{pert}}(\omega ;T) - \rho _{\mathrm{pert}}(\omega ;T=0), \end{aligned}$$
(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_{\mathrm{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 _{\mathrm{trans}} + \rho _{\mathrm{pert}}\), is employed as well. A variation of this Ansatz is also used in Ref. [22], 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,

$$\begin{aligned} K(\tau , \omega ) \rightarrow f(\omega ) K(\tau , \omega ), \qquad \rho (\omega ) \rightarrow \rho (\omega )/f(\omega ), \end{aligned}$$
(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 \rightarrow 0\) [16].

The Maximum Entropy Method (MEM) provides such a reconstruction method, based on Bayes’ theorem. The inversion is regularised by the prior probability (see below) and the most probable spectral function is obtained. In Bryan’s method [47], the (rescaled) spectral function is parametrised as [16, 48]

$$\begin{aligned} \frac{\rho (\omega )}{\omega } = \frac{m(\omega )}{\omega } \exp \sum _{i=1}^{N_{\mathrm{coeff}}} c_i u_i(\omega ), \end{aligned}$$
(28)

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

$$\begin{aligned} \int _0^\infty \frac{d\omega }{2\pi } \, u_i(\omega )u_j(\omega ) = \delta _{ij}. \end{aligned}$$
(29)

Since \(N_{\mathrm{coeff}} \sim N_\tau /2 \ll N_\omega \), this parametrisation leads to a significant reduction in the number of parameters to be determined. The form (28) is motivated by positivity, \(\rho (\omega )/\omega \ge 0\). In MEM, \(m(\omega )\) is referred to as the default model, see below. The conductivity is now determined by

$$\begin{aligned} \sigma \sim m^\prime (0) \exp \sum _i c_i u_i(0), \end{aligned}$$
(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)\) [48]. Here D indicates the data and H additional prior knowledge. The method relies on Bayes’ theorem,

$$\begin{aligned} P(\rho |DH) = \frac{P(D|\rho H) P(\rho |H)}{P(D|H)}, \end{aligned}$$
(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,

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

with

$$\begin{aligned} S(\rho ) = \int _0^{\infty } \frac{d\omega }{2\pi }\,\left[ \rho (\omega ) - m(\omega ) - \rho (\omega ) \ln \frac{\rho (\omega )}{m(\omega )} \right] , \end{aligned}$$
(33)

giving the method its name. The conditional probability now reads

$$\begin{aligned} P(\rho |DH) \propto e^{- L(\rho ) + \alpha S(\rho )}, \end{aligned}$$
(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. [48].

MEM has been applied to the conductivity in Refs. [15,16,17,18, 20, 21]. 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_{\mathrm{coeff}}\times N_\omega \) matrix, where \(N_{\mathrm{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. [20] and especially Ref.  [21] 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, in particular when the number of temporal points is sufficiently large. It is noted that the standard default model \(m(\omega )\sim \omega ^2\) for large \(\omega \) is modified due to perturbative corrections. On the lattice there are additional lattice artefacts due to the finite Brillouin zone [38, 46]. These effects can be incorporated in the reconstruction, by using more elaborate default models, see e.g. Ref. [43], and should lead to more sensitivity. Another 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 Bayesian methods, such as Bayesian Reconstruction (BR, see Sect. 3.7) [50, 51], or fully independent approaches, as discussed next.

3.5 Backus–Gilbert method

The Backus–Gilbert method is such an independent approach, designed for solving linear ill-defined problems with arguably controllable regularisation and systematic uncertainty, although this may be hard to quantify in practice. Rather than reconstructing the entire spectral function \(\rho (\omega )\), it aims to represent it by an estimator

$$\begin{aligned} {{\widehat{\rho }}}(\omega _0) = \int _0^\infty d\omega \, \delta (\omega _0, \omega ) \rho (\omega ), \end{aligned}$$
(35)

where \(\delta (\omega _0, \omega )\) is called the resolution function, which should be narrowly peaked around \(\omega _0\) and normalised,

$$\begin{aligned} \int _0^\infty d\omega \,\delta (\omega _0, \omega ) = 1, \end{aligned}$$
(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 [22, 24]

$$\begin{aligned} \delta (\omega _0, \omega ) = \sum _{n=1}^{N_\tau - 1} q_n(\omega _0) K(\tau _n, \omega ). \end{aligned}$$
(37)

The functions \(q_n(\omega _0)\) are found by minimising the second moment of the resolution function squared,

$$\begin{aligned} \varGamma _{\omega _0} = \int _0^\infty d\omega \, (\omega - \omega _0)^2 \delta ^2(\omega _0, \omega ), \end{aligned}$$
(38)

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

$$\begin{aligned} q_n(\omega _0) = \frac{W(\omega _0)^{-1}_{nm} R_m}{R_k W(\omega _0)^{-1}_{kl} R_l}, \end{aligned}$$
(39)

in terms of

$$\begin{aligned} W(\omega _0)_{nm}= & {} \int _0^\infty d \omega \, (\omega - \omega _0)^2 K(\tau _n, \omega ) K(\tau _m, \omega ), \end{aligned}$$
(40)
$$\begin{aligned} R_n= & {} \int _0^\infty d \omega \, K(\tau _n, \omega )\,. \end{aligned}$$
(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

$$\begin{aligned} {{\widehat{\rho }}}(\omega _0)= & {} \int _0^\infty d\omega \, \sum _n q_n(\omega _0) K(\tau _n,\omega )\rho (\omega ) \nonumber \\= & {} 2\pi \sum _n q_n(\omega _0) G(\tau _n). \end{aligned}$$
(42)

The conductivity is then extracted in a usual way,

$$\begin{aligned} \sigma \sim \lim _{\omega _0 \rightarrow 0} \frac{{{\widehat{\rho }}}(\omega _0)}{\omega _0}, \end{aligned}$$
(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

$$\begin{aligned} {{\widehat{\rho }}}(\omega _0) = f(\omega _0) \int _0^\infty d \omega \, \delta (\omega _0, \omega ) \frac{\rho (\omega )}{f(\omega )}, \end{aligned}$$
(44)

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

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

where \(S_{nm}\) is the covariance matrix for the correlator \(G(\tau _n)\), and \(\lambda \) is a tunable parameter, determined empirically 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. [22, 24].

3.6 Tikhonov regularisation

The method of additive regularisation, see Eq. (45), is not the only one. Tikhonov regularisation [49] 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

$$\begin{aligned} W= & {} U \varSigma V^T,\quad U^T U = V^T V = 1\!\!1, \nonumber \\ \varSigma= & {} \text{ diag }(\sigma _1,\,\sigma _2,\,\ldots ,\,\sigma _N), \end{aligned}$$
(46)

with \(\sigma _1 \ge \sigma _2 \ge \ldots \ge \sigma _N\). Since \(W^{-1} = V \varSigma ^{-1} U^T\), the inversion of W comes down to the inversion of \(\varSigma \), which can be regularised as follows,

$$\begin{aligned} \varSigma ^{-1} = \text{ 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) , \end{aligned}$$
(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. [24].

Table 1 Details of the lattice QCD ensembles to compute the electrical conductivity. Fermion properties refer to sea quarks. Here \(a_\tau \) and \(a_s\) denote the temporal and spatial lattice spacing respectively
Full size table

3.7 Other approaches

Here we briefly list some additional inversion methods. References  [28, 29] further develop the proposal of Ref. [27], 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 Sect. 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 Reconstruction (BR) approaches in an extended search space can be found in Refs.  [50, 51]. Ref.  [52] 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. [52] 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. [53]. In Ref.  [54] kernel ridge regression was applied to mock data in quantum many-body physics; a first application to QCD data can be found in Ref. [55] 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 Sect. 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. References [15,16,17,18, 23] concern quenched QCD, while in Refs. [19,20,21,22, 24] 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.  [24], 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. [20, 21], 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 \rightarrow 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 2 Lattice sizes used compute the electrical conductivity: number of spatial (\(N_s\)) and temporal (\(N_\tau \)) lattice points, and corresponding temperatures
Full size table
Table 3 Details of the current and inversion method used to compute the electrical conductivity
Full size table

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.  [20, 21], 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 [16, 24],

    $$\begin{aligned} G_{ij}(x)= & {} \langle A_i(x) A_j(0) \rangle - (-1)^{\tau /a_\tau } \langle B_i(x) B_j(0) \rangle , \nonumber \\ A_i= & {} {{\bar{\psi }}}\gamma _i \psi , \qquad \quad B_i = {{\bar{\psi }}} \gamma _4 \gamma _5 \gamma _i \psi . \end{aligned}$$
    (48)

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

    $$\begin{aligned} \rho _{ij}(\omega ) = \frac{1}{2}\left[ \rho ^{\mathrm{even}}_{ij}(\omega ) + \rho ^{\mathrm{odd}}_{ij}(\omega ) \right] , \end{aligned}$$
    (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^{\mathrm{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 [15, 16] 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)

    $$\begin{aligned} j^{\mathrm{cons}}_{x,i}&= \frac{1}{2}{{\bar{\psi }}}_{x+{{\hat{\imath }}}}(1+\gamma _i) U^\dagger _{x,i}\psi _x\nonumber \\&- \frac{1}{2}{{\bar{\psi }}}_x(1-\gamma _i) U_{x,i}\psi _{x+{{\hat{\imath }}}}, \end{aligned}$$
    (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.  [20, 21, 24], 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. [22] employs a mixed combination, \(\langle j^{\mathrm{cons}}_{x,i}j^{\mathrm{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^{\mathrm{loc}}_{x,i}j^{\mathrm{loc}}_{y,i}\rangle /\langle j^{\mathrm{cons}}_{x,k}j^{\mathrm{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 Sect. 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 \({{{\mathcal {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. Refs. [56, 57]) and can provide an additional handle on hydrodynamic behaviour at small \(\omega \) and \(|{{\mathbf {p}}}|\), as well as a better constraint of the spectral function by using perturbative information.

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_{\mathrm{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.

Fig. 2
figure 2

Results for the electrical conductivity, normalised as \(\sigma /(TC_{\mathrm{em}})\), in quenched QCD (\(N_f=0\)) as a function of \(T/T_c\). The early result [15], \(\sigma /(T C_{\mathrm{em}}) \approx 7\), is not shown for clarity

Full size image

Results in quenched QCD are shown in Fig. 2, as a function of \(T/T_c\). The very early result from Ref. [15], \(\sigma /(T C_{\mathrm{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 /(T C_{\mathrm{em}}) \approx 0.2-0.5\). Note that Ref. [17] provides both a precise result, \(\sigma /(TC_{\mathrm{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. [16] on the one hand and Refs. [17, 18, 23] on the other hand, making the agreement is perhaps less surprising. In any case, it is interesting that the early quenched result of Ref.  [16], obtained using staggered quarks without taking a continuum limit, remains to be consistent with the renormalised continuum-extrapolated Wilson-clover results of Refs.  [17, 18, 23]. 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.

We note here that the error bars are dominated by uncertainties due to the reconstruction method. Statistical uncertainties in the current-current correlators themselves are under control. To assess the uncertainty in the reconstructed spectral functions, the statistical uncertainty is typically evaluated using a jackknife or bootstrap analysis. Estimates of systematic uncertainties are method-dependent. These include using different model functions (for fits), and varying default models (for MEM) and regularisation parameter \(\lambda \) (for Backus-Gilbert). In all methods, the rescaling function \(f(\omega )\) – see Eq. (27) – and the number of time slices included can be varied to investigate robustness. With regard to the summary of results in Figs. 2 and 3, it is important to keep in mind that each publication partly follows its own method-dependent strategy to estimate the final error, which can differ quite significantly; we refer to the original publications for further detail and come back to this in the Conclusion.

Fig. 3
figure 3

Results for electrical conductivity, normalised as \(\sigma /(TC_{\mathrm{em}})\), in QCD with \(N_f=2\) and \(2+1\) dynamical flavours as a function of temperature in MeV

Full size image

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. [24] follows the lattice formulation and choice of parameters of Refs. [58, 59]. Moreover, Ref. [24] 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. [20,21,22] are higher than in Ref. [24], 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. [20, 21] and Ref. [22], 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. [20,21,22], compared to Ref. [24], 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. [19, 22] 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_{\mathrm{em}})\) is comparable to the one obtained in the quenched case.

Table 4 Estimates of the pseudocritical temperatures for the studies with dynamical quarks
Full size table

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_{\mathrm{em}})\sim 1/g^4\ln 1/g\) [9], 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 \({{{\mathcal {N}}}} = 4\) Yang- Mills theory at nonzero temperature [13, 14]. In Ref. [21] 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. [21].

A related aspect is the width of the so-called transport peak, assuming it is present. As mentioned in Sect. 2, in the case of a narrow transport peak, with a width \(\gamma \ll T\), indicative of weak coupling, the correlator is only sensitive to its area but not to further details, making the extraction of the conductivity not possible. Concerning the work discussed above, it is possible to analyse the observed widths, either directly from the spectral functions reconstructed or because the width is a fit parameter. It is observed that the transport peak is indeed not narrow, and the widths fall in the range \(\gamma \sim 1-4T\). In this case the correlator is sensitive to details of the peak, and assuming, e.g. in the case of a fit, that the peak is still characterised by its height and width, the conductivity is accessible. A caveat here is that this reasoning is at best self-consistent (or perhaps self-fulfilling), since only broad transport peaks lead to accessible conductivities. The insensitivity to narrow peaks [39] remains an outstanding open problem.

Before concluding this section, we note that the lattice data of Ref. [17] have been re-analysed in two papers. Ref. [29] employed the approach of Ref.  [27], see also Ref. [28], 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.  [45], 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. [17], namely, at \(T/T_c=1.45\),

$$\begin{aligned} \sigma /(T C_{\mathrm{em}}) > rsim ~ 0.1, \qquad \quad 2\pi T D > rsim 0.8, \end{aligned}$$
(51)

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

The data of Refs. [17] and [22] has also been re-analysed in Refs. [25] and [26] respectively, using thermal sum rules to constrain the Ansätze used in the fits. In Ref.  [25], a higher result was found for quenched study at \(T/T_c=1.45\), namely

$$\begin{aligned} \sigma /(T C_{\mathrm{em}}) \sim 0.57. \end{aligned}$$
(52)

In Ref. [26] the \(N_f=2\) data given in Ref.  [22] 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.  [60, 61], 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. [62] 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

$$\begin{aligned} \frac{\sigma (\mu )}{\sigma (0)} = 1+ c(T)\frac{\mu ^2}{T^2} + {{{\mathcal {O}}}}\left( \frac{\mu ^4}{T^4}\right) , \end{aligned}$$
(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 [24]. 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 [5] 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. One way forward is to contrast and combine the various approaches for spectral reconstruction, using the same data. In addition, progress can be made by developing a common framework for estimates of systematic uncertainties in spectral reconstruction, especially since each approach comes with its own methodology, making a direct comparison more involved. It is also noted that none of the results include the disconnected contributions yet, for reasons discussed in Sect. 4. It would be worthwhile to estimate the importance of those eventually.

Nevertheless, a comparison between the existing lattice studies, presented in Sect. 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 /(T C_{\mathrm{em}}) \lesssim 0.5\), where \(C_{\mathrm{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 /(T C_{\mathrm{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.  [63] 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 [64]. Another interesting possibility is the analysis of the current-current correlator at imaginary chemical potential.