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

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\"atze 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.


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 Nuclotronbased 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 σ 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 a Email: g.aarts@swan.ac.uk b Email: aleksandr.nikolaev@swan.ac.uk 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 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 τ 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 [25][26][27] and neither other transport coefficients such as the shear [28][29][30] and bulk viscosities [31,32]. 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 σ 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.

Kubo formula and spectral function
The electromagnetic current in QCD receives contributions from all quark flavours and reads 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 em µ † (x) = j em µ (x). The electrical conductivity σ indicates the linear relationship between the current density and an electric field, j em i = σ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, at vanishing energy and momentum, i.e., 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, where χ Q is the charge susceptibility, Here V is the spatial volume and Q is the total charge, i.e. the volume integral of j em 0 (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, via the standard relation [see Eq. (69)] with the kernel Here the Euclidean time 0 ≤ τ < 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 p dependence. When prefactors involving eq f are dropped, the superscript 'em' is omitted as well.
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 [33] with N c = 3. Here n F (ω) = 1/[exp(ω/T )+1] is the Fermi-Dirac distribution. The quantity I in the first term reads with ω k = √ k 2 + m 2 . For massless quarks this evaluates as 1 Note there is a typo in the last line of Eq. (19) in Ref. [33]: H should read a H .
the spectral function is odd, ρ ii (−ω) = −ρ ii (ω), as it should be. Below we take ω ≥ 0. The corresponding Euclidean correlator reads, in the massless limit [33], where u = 2πT (τ −1/2T ). The first (constant) term comes from the first term in Eq. (9); the τ dependent term from the second one. At the midpoint, τ = 1/2T , the contributions from both terms are comparable, 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 ω = ω k + ω p+k (with k the loop momentum and p = 0), which is permissible once |ω| > |2m|. It contains the vacuum contribution, increasing as ω 2 at large ω, 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 ω 3 as ω → 0. We will refer to this term as the continuum or perturbative contribution.
The first term, with ωδ(ω), is only present at nonzero temperature and arises from the cut corresponding to scattering with particles in the heatbath, ω k + ω = ω p+k , in the limit that p → 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 δ function in Eq. (9) is smeared out and takes the form of a so-called transport peak, where γ is proportional to the collisional scattering width or the inverse mean free path, A trans (= 2N c I) is the overall coefficient, and σ ∝ A trans /γ. The origin of the transport peak can be seen in various ways, using e.g. Feynman diagrams and pinching poles [34,35] or kinetic theory [36,37]. 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 (γ 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 [34]. In the limit that ω < Λ T , the kernel (8) simplifies to 2T /ω, and integrating the transport peak (14) in Eq. (7) yields where in the last expression the cutoff Λ on ω has been removed. The crucial observation is that this expression is independent of γ and the Euclidean time τ , indicating the insensitivity of the correlator to narrow transport peaks. In fact, taking A 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 =ψγ i ψ, with itself. In lattice QCD, this is commonly referred to as the disconnected contribution. Perturbatively, they start contributing at O(α 3 s ) 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 In the disconnected contribution on the other hand, we find the square of the sums, 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.
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.

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 τ N τ . Here a τ is the temporal lattice spacing and N τ the number of points in the time direction; Euclidean time is discretised as τ /a τ = 0, . . . , N τ − 1. Since the correlator G(τ ) is known at a finite number of temporal points and the spectral function ρ(ω) is in principle a continuous function of ω, this inversion problem is far from straightforward. To be more precise, due to reflection symmetry, K(τ, ω) = K(1/T − τ, ω), the number of points available for the analysis is on the order of N τ /2; even after placing an upper limit on the ω interval, such that 0 < ω < ω max , and discretising the finite interval, typically on the order of N ω = 1000 points are used to present ρ(ω). Since N ω N τ , the inversion problem is ill-posed. In addition, the focus on the ω → 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.

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(τ, ω), due to the compact time direction, 0 ≤ τ < 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 [38,39]. Let us suppose a spectral function ρ(ω; 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 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).

Sum rules
Exact sum rules are important [19] to constrain the currentcurrent spectral function at nonzero temperature. Defining the difference between the finite and zero-temperature spectral function as one finds the sum rule [19], in the thermodynamic limit, Note that the zero-temperature ω 2 contribution cancels in the subtraction (19). Since the Operator Product Expansion (OPE) predicts [40] that the thermal contribution decays as (T /ω) 2 at large ω, 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. [41,42].

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 ρ(ω) 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 ω, with in particular a linear slope in ω; -continuum (ω 2 ) contribution at high ω, possibly modified by lattice artefacts [33,43]; -at least one bound-state peak in the low-temperature phase, to represent the vector meson.
Refs. [17,18,23] employ an Ansatz combining a transport peak and the expected perturbative continuum behaviour (for massless quarks) in the deconfined phase, where ρ trans (ω) is given in Eq. (14) and c.f. Eq. (9). The three temperature-dependent parameters are the coefficients A trans,pert and the width γ of the transport peak. Note that A 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 σ [23].
Ref. [19] employs a related Ansatz for the subtracted spectral function, ∆ρ(ω, 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 with the new term Here m B , g B and A bound indicate the mass, width and strength of the bound state. The factor tanh(ω/T ) 3 ensures the contribution does not contribute to the conductivity in the ω → 0 limit and decays as 1/ω 2 at large ω, as predicted by the OPE [40]. It is also noted in Refs. [19,44] that the transport peak (14) in fact violates this condition.
Hence it is proposed [19] to modify it as which still has linear behaviour at small ω but decays as 1/ω 2 at large ω. Finally, the main reason for the subtraction is to eliminate the zero-temperature ω 2 contribution, 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 ω 2 term. Overall, the number of parameters (A trans,pert,bound , m B , g B , γ) 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 ρ trans + ρ 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.

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 illposedness of the inversion. Before proceeding, we note the possibility to rescale the kernel and the spectral function, (27) leaving the product unchanged, to stabilise the inversion. A common rescaling is to use f (ω) = ω, to resolve the 1/ω divergence in the kernel as ω → 0 [16].
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 [45], the (rescaled) spectral function is parametrised as [16,46] with u i (ω) (i = 1, . . . , N coeff ) an orthogonal but incomplete set of basis functions, The reduction follows since 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 (ρ|DH) [46]. Here D indicates the data and H additional prior knowledge. The method relies on Bayes' theorem, where P (A|B) stands for the conditional probability of A given B. In this expression, P (D|ρH) is the likelihood function, P (ρ|H) the prior probability, and P (D|H) a normalisation. While the likelihood function, P (D|ρH) = e −L(ρ) , is familiar from standard χ 2 -minimisation, the prior probability contains an entropy-like term, with giving the method its name. The conditional probability now reads P (ρ|DH) ∝ e −L(ρ)+αS(ρ) , with α determining the balance between the two terms. While at α = 0 the method reduces to a standard fitting procedure, in absence of any data the probability is extremised when ρ(ω) = m(ω), yielding the default result. For further details on MEM we refer to Ref. [46]. MEM has been applied to the conductivity in Refs. [15][16][17][18]20,21]. We note here that the basis functions u i (ω) are obtained via a singular value decomposition (SVD) of the kernel K(τ, ω), when viewed as a N coeff ×N ω matrix, where N coeff N τ /2, linking the size of the set of basis functions to the temporal extent. This is a limitation for small N τ 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(ω) is quite mild [21]. 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.

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 ρ(ω), it aims to represent it by an estimator where δ(ω 0 , ω) is called the resolution function, which should be narrowly peaked around ω 0 and normalised, 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] δ(ω 0 , ω) = The functions q n (ω 0 ) are found by minimising the second moment of the resolution function squared, which should effectively minimise its width. Combining the equations above, one obtains the solution (summation over repeated indices implied) in terms of Note that since the kernel K(τ, ω) is symmetric around the midpoint, the number of independent basis functions q n is limited by N τ /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 The conductivity is then extracted in a usual way, assuming that the estimator ρ(ω 0 ∼ 0) is close to the physical spectral function. The goal is therefore to find the optimal set of functions q n (ω 0 ) which make δ(ω 0 , ω) as narrow as possible in ω for ω 0 ∼ 0. The main difficulty in this formulation is the inversion of the matrix W (ω 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 and finding an optimal choice for f (ω) (e.g. f (ω) = ω). Secondly, the matrix W can be regularised as where S nm is the covariance matrix for the correlator G(τ n ), and λ is a tunable parameter, determined by comparing the behaviour of δ(ω 0 , ω) for different values of λ. Further discussion on successes and limitations of this approach in the context of the conductivity can be found in Refs. [22,24].

Tikhonov regularisation
The method of additive regularisation, see Eq. (45) Σ = diag(σ 1 , σ 2 , . . . , σ N ), with σ 1 ≥ σ 2 ≥ . . . ≥ σ N . Since W −1 = V Σ −1 U T , the inversion of W comes down to the inversion of Σ, which can be regularised as follows, where the parameter has to be chosen carefully; small 's lead to precise but unstable results, while large 's guarantee stable inversion at a cost of loss of accuracy. Further experiments with this approach can be found in Ref. [24].

Other approaches
Here we briefly list some additional inversion methods. Refs. [44,49] further develop the proposal of Ref. [48], 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 shortdistance 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. [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.

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. [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 τ (the spatial and temporal lattice spacing respectively), except Refs. [20,21], in which anisotropic lattices, with a s /a τ = 3.5, are employed. The advantage of the latter is the finer temperature resolution, when N τ is varied. The lattice cutoff is dealt with in a variety of ways: continuum limit, indicated by a τ → 0; -fixed scale: the lattice cutoff is fixed and temperature is varied by changing N τ , according to T = 1/(a τ N τ ); -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 τ < 0.0618 fm. Table 2 contains the details of the lattice geometry, i.e. the number of lattice points in spatial (N s ) and temporal (N τ ) 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], To resolve this staggering, the spectral reconstruction is performed on even and on odd time slices independently, obtaining two spectral functions, ρ even,odd ij (ω). In the sum the oscillating contribution cancels, and one can proceed with ρ ij (ω) as usual. This procedure limits, however, the number of time slices effectively available by a factor of two. This issue does not arise with Wilsontype fermions, which are hence the preferred choice. -local/conserved current: the simplest choice for the current operator is the local one, j loc x,i =ψ x γ i ψ 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.  [17] 1012.4963 0 quenched − aτ → 0 1 continuum limit [18] 1112.4802 0 quenched − 0.015 1 fixed scale [19] 1212.4200 2 Wilson-clover 270 0.0486(4)(5) 1 fixed cutoff [20] 1307.6763 2 + 1 Wilson-clover 384(4) 0.0350 (2) 3.5 fixed scale [21] 1412.6411 2 + 1 Wilson-clover 384(4) 0.0350 (2) 3.5 fixed scale [22] 1512.07249 2 Wilson-clover 270 0.0486(4)(5) 1 fixed scale [23] 1604.06712 0 quenched − aτ → 0 1 continuum limit [24] 1910.08516 2 + 1 staggered 134.2(6) 0.0618, 0.0493 1 fixed cutoff  Table 3. Details of the current and inversion method used to compute the electrical conductivity. ref.
fermion type current renormalised inversion method [15] staggered local − Bayesian priors, MEM, Ansatz [16] staggered local − MEM [17] Wilson-clover local Ansatz, MEM [18] Wilson-clover local Ansatz, MEM [19] Wilson-clover local Ansatz, sum rule constraints [20,21] Wilson-clover conserved MEM [22] Wilson-clover mixed local-conserved Ansatz, sum rule constraints, Backus-Gilbert [23] Wilson-clover local Ansatz [24] staggered conserved Backus-Gilbert, Tikhonov regularisation 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) 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, j x,i j y,i , 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, j cons x,i j loc y,i , 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 j loc x,i j loc y,i / j cons x,k j loc y,k . -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 O(α 3 s ) 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. [56]) and can provide an additional handle on hydrodynamic behaviour at small ω and |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 em = 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.  Results in quenched QCD are shown in Fig. 2, as a function of T /T c . The very early result from Ref. [15], σ/(T C em ) ≈ 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 σ/(T C em ) ≈ 0.2 − 0.5. Note that Ref. [17] provides both a precise result, σ/(T C 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 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. [57,58]. Moreover, Ref. [24] has results at two lattice spacings; the red star symbols denote the results at the smaller spacing, N τ = 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 σ/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 σ/(T C 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 σ/(T C em ) ∼ 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 = σ/χ Q , where χ Q is the charge susceptibility. The characteristic result at strong coupling is D = 1/(2πT ) in 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 χ Q within the same lattice QCD setup, with the result that 0.5 < 2πT D < 2, compatible with the holographic order of magnitude at strong coupling. Moreover, it was observed that 2πT D has a minimum in the crossover region, see Fig. 14 of Ref. [21].
Before concluding this section, we note that the lattice data of Ref. [17] have been re-analysed in two papers. Ref. [44] employed the approach of Ref. [48], see also Ref. [49], 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. [40], 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, which is indeed smaller by a factor of 3 for the conductivity. It is stated [44] 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. [41] and [42] respectively, using thermal sum rules to constrain the Ansätze used in the fits. In Ref. [41], a higher result was found for quenched study at T /T c = 1.45, namely σ/(T C em ) ∼ 0.57.
In Ref. [42] 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.

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. [59,60], 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. [61] 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 µ, the dependence on chemical potential is considered via the expansion and it is found that the maximal value of the second-order coefficient, c(T ) ≈ 0.15 (5), is reached in the vicinity of the chiral crossover. Hence the coefficient c(T ) is quite small, and even at µ/T ≈ 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 µ 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 B = 0 results of this reference have been discussed above; with nonzero 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.

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 σ/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 σ/(T C em ) 0.5, where C 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 σ/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 σ/(T C 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 = σ/χ Q , where χ 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. [62] 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 µ/T is possible, although it is expected to be noisy numerically [63]. Another interesting possibility is the analysis of the currentcurrent correlator at imaginary chemical potential.

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 , x). These relations are well known [64,65].
Let us start with the retarded and advanced Green's functions G R (x−y) = iθ(t x −t y ) [O(x), O † (y)] = G A (y −x), (54) and the spectral function Expectation values are taken in thermal equilibrium, which explains the x − y dependence. After going to momentum space and using the identity we arrive at the dispersion relation Employing the identity then yields the important relation, ρ(ω, p) = −i [G R (ω, p) − G A (ω, p)] = 2Im G R (ω, 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, with 0 ≤ τ < 1/T , is written in momentum space as where ω n = 2πnT , n ∈ Z, are the Matsubara frequencies (we consider bosonic operators here). By analyticity, it satisfies a similar dispersion relation as above, leading to the important relation If a Euclidean correlator is known analytically, the spectral function can be obtained following the sequence G E (τ, x) → G E (ω n , p) → G R (ω, p) → ρ(ω, p).
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), where we used the identity n B (ω) + n B (−ω) + 1 = 0.