Electrical conductivity and charge diffusion in thermal QCD from the lattice

We present a lattice QCD calculation of the charge diffusion coefficient, the electrical conductivity and various susceptibilities of conserved charges, for a range of temperatures below and above the deconfinement crossover. The calculations include the contributions from up, down and strange quarks. We find that the diffusion coefficient is of the order of 1/(2\pi T) and has a dip around the crossover temperature. Our results are obtained with lattice simulations containing 2+1 dynamical flavours on anisotropic lattices. The Maximum Entropy Method is used to construct spectral functions from correlators of the conserved vector current.


Introduction
Recently, there has been a great deal of progress in the understanding of the dynamical and static properties of the quark-gluon plasma (QGP). The current theoretical interpretation of heavy-ion collision experiments at RHIC and CERN consists of a hydrodynamical description of the evolution of the fireball, see e.g. the reviews [1,2]. This relies on the QGP medium thermalising after a very short time of less than 1 fm/c; the subsequent evolution is then modelled by viscous hydrodynamics, until hadronisation.
The input parameters in the hydrodynamic evolution equations are the equation of state and, for the nonideal case, transport coefficients, such as viscosities and conductivities. These quantities capture the dynamics from the underlying theory and hence are determined by QCD. A first-principles determination must face the challenge of strong coupling: it is now widely accepted that dynamics in the QGP is strongly coupled, typically expressed via the statement that the ratio of the shear viscosity to entropy density η/s is close to the value obtained in models with holographic duals [3].
Besides the shear and bulk viscosities, there is an increasing interest in another transport coefficient, namely the electrical conductivity, due to the role it plays in heavy-ion phenomenology. For instance, the conductivity has recently been discussed in the context of charge density fluctuations [4] and the evolution of strong electromagnetic fields produced in noncentral collisions [5][6][7][8][9][10]. It has also been suggested that experimental information on conductivity can be extracted from flow parameters in heavy-ion collisions [11].
Using linear response theory, transport coefficients can be related to current-current spectral functions in thermal equilibrium, via the celebrated Kubo relations [12,13]. This opens up the possibility to apply lattice QCD at finite temperature as the nonperturbative tool of choice, provided that the analytical continuation from Euclidean to real time, or from Euclidean correlators to spectral functions, can be carried out reliably [14]. In recent years several results have been obtained along these lines. Refs. [15,16] contain the first (theoretically controlled) results for the shear and bulk viscosity, see also the review [17]. The best-studied transport coefficient is however the electrical conductivity, since the corresponding Euclidean correlator can be computed numerically with high precision. The first results were obtained using the staggered fermion formulation [18,19] in quenched QCD. Since then the study of the systematic uncertainties and the extension to Wilsontype quarks have taken a central role, in quenched QCD [20,21] and in dynamical QCD with N f = 2 and 2 + 1 flavours [22,23]. All studies [19][20][21][22][23] are in qualitative agreement: around T = 1.5T c , where T c is the crossover temperature, the conductivity is of the order of σ = 0.2 − 0.4C em T c , where C em is an electromagnetic prefactor depending on the quark charges (see Sec. 5 below). Recent non-lattice studies include Refs. [24][25][26][27][28][29][30][31].
The most detailed lattice study so far can be found in Ref. [23], where the temperature dependence of σ/T across the deconfinement transition was studied for the first time, over a range of temperatures corresponding to 0.63 − 1.9T c . In that analysis we used a QGP with N f = 2 + 1 flavours but only included the light quark contribution to the conserved vector current. Here we improve upon those results by also including the strange quark contribution and comparing the relative contributions. Moreover we compute various susceptibilities, including the charge susceptibility χ Q , which allows us to compute for the first time the charge diffusion coefficient D = σ/χ Q in a self-contained calculation. Since the diffusion coefficient can also be computed in strongly coupled theories that permit a holographic prescription, with the characteristic result that D = 1/(2πT ) in N = 4 Yang-Mills theory at nonzero temperature [32][33][34], this direct computation allows us to compare QCD with strongly coupled gauge theories which have a dual formulation.
The remainder of the paper is organised as follows. In the next section we start with the information on the lattice action and ensembles used in this work, followed by a determination of the crossover transition temperature using the renormalised Polyakov loop and the chiral susceptibility in Sec. 3 and the baryon, isospin and charge susceptibilities in Sec. 4. In Sec. 5 we turn to the electrical conductivity, and present the Euclidean correlators and their corresponding spectral functions determined via the Maximum Entropy Method. The systematics in this construction are discussed in some detail. Finally, our results for the diffusion coefficient, obtained using the results from the two preceding sections, are presented in Sec. 6. We summarise our findings and provide a brief outlook in Sec. 7.

N f = + 1 lattice details
In this work we follow the Hadron Spectrum Collaboration [35] and use a Symanzikimproved anisotropic gauge action with tree-level mean-field coefficients and a mean-fieldimproved Wilson-clover fermion action with stout-smeared links [36]. The anisotropy, with a reduced temporal lattice spacing a τ < a s , is crucial to obtain a better resolution of the correlation functions, especially at higher temperatures. This will be discussed further below. Anisotropy introduces two new bare parameters in the action, the bare gauge and fermion anisotropies, which are nonperturbatively tuned to give the desired renormalised anisotropy, ξ = a s /a τ , common to the gauge and fermionic degrees of freedom. The ensembles employed here are part of our "2nd generation" data set [37] and were previously used for a determination of the conductivity (from two light flavours only) [23] and the bottomonium spectrum at nonzero temperature [38]. The gauge action takes the form (2.1) where P and R are the 1 × 1 plaquette and 2 × 1 rectangular Wilson loops, u s(τ ) are the spatial (temporal) mean links, γ g is the bare gauge anisotropy and, as usual, β = 2N c /g 2 with N c = 3 colours.
The fermion action (for a single flavour) reads wherem 0 = a τ m 0 is the bare fermion mass, γ f the bare fermion anisotropy, ∇ µ covariant finite differences, σ µν = i 2 [γ µ , γ ν ], and the clover coefficients The spatial gauge links in the fermion action have been stout smeared [36] with smearing weight ρ = 0.14 and n ρ = 2 iterations, andũ s is the mean value of the spatial stout-smeared links.
The choice of bare parameters is given in Table 1 and follows the tuning by the Hadron Spectrum Collaboration [35]. The resulting renormalised anisotropy is ξ = 3.5. The two  Table 2. Details of the ensembles. The lattice size is N 3 s ×N τ , with the temperature T = 1/(a τ N τ ). N CFG (N SRC ) denote the number of configurations at each volume (the number of source positions within the volume) used for the analysis of the conductivity. degenerate light quarks yield M π = 384(4) MeV (2.8 times larger than the physical pion), corresponding to M π /M ρ = 0.446(3) [39], while the third flavour is tuned to the strange quark mass [35].
Details of the finite-temperature ensembles are given in Table 2. Note that there are five ensembles in the hadronic phase and four in the quark-gluon plasma phase. The determination of the pseudo-critical temperature T c is discussed in the next section. In order to look for finite-size effects, we have generated configurations with two spatial volumes at four different temperatures, with spatial extents of ∼ 2.9 respectively 3.9 fm (N s = 24 and 32). Access to the zero-temperature configurations (N τ = 128) has been kindly provided to us by the Hadron Spectrum Collaboration. The ensembles were generated using the Rational Hybrid Monte Carlo (RHMC) algorithm with multiple timescale integration and Hasenbusch preconditioning for the light quarks, using the Chroma software suite [40] with Bagel routines [41]. For further details about the algorithm, we refer to Refs. [35,39]. After 1000 thermalisation trajectories (2000 for the 32 3 × 24 ensemble), configurations were sampled every 10 RHMC trajectories, except for the 32 3 × 48 ensemble where the sampling frequency was every 5 trajectories. The plaquette and Polyakov loop autocorrelation times were found to be between 2 and 30 trajectories.

Deconfinement and chiral transition
After a determination of the lattice spacing, the temperature can be specified in MeV very precisely using the standard relation T = 1/(a τ N τ ). However, it is desirable to express the temperature in units of T c , the crossover temperature, especially since the two light quarks are heavier than in Nature (for a lattice study of the transition with physical quark masses, see e.g. Ref. [42]). In order to do so, we use the renormalised Polyakov loop as an indicator of the deconfinement transition, following closely the renormalisation procedure described in Ref. [43].
The Polyakov loop expectation value L is related to the free energy of a static quark F via However, F is only defined up to an additive renormalisation constant ∆F , which depends on the gauge coupling and other bare parameters but not on the temperature. Expressing the renormalised free energy F R as F R = F +∆F allows us write the renormalised Polyakov loop as which defines the multiplicative renormalisation constant Z L . Following Ref. [43], we impose a renormalisation condition at a reference temperature T R , by requiring that which fixes Z L . Figure 1 shows the Polyakov loop with three different renormalisation schemes corresponding to different choices of T R and the constant in Eq. (3.3), as detailed in the figure caption. The data are interpolated using cubic splines, with the statistical uncertainty given by the thickness of the three interpolating curves; it can be seen to be negligible. We then obtain the Polyakov loop susceptibility as the derivatives of the interpolating curves for the Polyakov loop in each of the three schemes. The peak positions are indicated with the vertical lines in Fig. 1 and give us the point of inflection at N crit (7), where the error reflects the systematic error coming from the spread of the three renormalisation schemes. This corresponds to a deconfinement critical temperature of T c = 185(4) MeV. We note that neither chiral nor continuum extrapolations have been performed in our analysis.
In the limit of massless quarks, QCD becomes classically invariant under chiral transformations. This symmetry is spontaneously broken at low temperature. For physical masses, even if the chiral symmetry is explicitly broken, the associated order parameter shows a clear transition signal at a certain temperature. The chiral condensate and, in particular, its susceptibility χ c are commonly used to define the crossover transition temperature. We have determined the chiral susceptibility due to the two light flavours, using [44][45][46] where Z is the partition function, M the fermion matrix, V the spatial volume and m the degenerate light quark mass. Moreover, we introduce here the connected χ conn and the disconnected χ disc contributions to the susceptibility. The traces in Eq. one. Because we change the temperature by changing the value of N τ rather than the lattice spacing, the (additive and multiplicative) renormalisation of χ c is the same for all temperatures. A peak in the susceptibility occurs therefore at the same temperature for the renormalised and unrenormalised χ c . We hence show the bare susceptibility and are only interested in the overall shape, rather than the absolute value. The results are plotted in Fig. 2. The connected contribution is not singular while the disconnected contribution shows a peak. Because the determination of the disconnected term is particularly expensive at low temperature, we could not determine the chiral critical temperature very precisely. Our best estimate is T χ c ≈ 170(20) MeV, which is somewhat lower than the value obtained from the Polyakov loop. In the following we will use the value for T c determined from the Polyakov loop as this has been obtained with a better precision.

Susceptibilities
Fluctuations of conserved charges are sensitive probes both of the thermal state of the medium and of its critical behaviour. They are quantified by susceptibilities, defined as second (and higher) derivatives of the free energy with respect to the chemical potential associated with the investigated charge. In QCD, assuming three active light flavours, charges that can be studied include baryon number, electric charge and strangeness. Their susceptibilities probe the actual degrees of freedom that carry such charges, i.e. quarks or hadrons. Experimentally, fluctuations can be used to probe quark confinement [47] by studying event-by-event fluctuations of charged particle ratios [48]. Susceptibilities show a rapid rise in the crossover region: at low temperature they are small since quarks are confined; at high temperature they are larger and they approach the ideal gas limit. They have been studied by many groups in the past [49][50][51][52][53][54]. Notably, so far, lattice studies have mainly employed staggered fermions. Here instead we use clover-improved Wilson fermions.
For an earlier study using Wilson fermions see Ref. [55]. The charge diffusion coefficient D and the electrical conductivity σ are related via the well-known relation D = σ/χ Q , where χ Q is the charge susceptibility [13]. In this section we determine χ Q and various other (second-order) susceptibilities, defined as second derivatives of the free energy with respect to the chemical potential associated with a conserved charge. We introduce the quark number density and the quark number susceptibilities where Z is the partition function, V the spatial volume, and µ f the quark chemical potentials for flavours f ∈ {u, d, s}. Baryon (B), isospin (I) and electrical charge (Q) chemical potentials are related to the quark chemical potentials as In general we denote the electrical charge of the quark as eq f , with e the elementary charge and q f = 2/3 or −1/3 its fractional charge. All desired quantities can now be expressed in terms of n f and χ f f . The baryon number density and baryon number susceptibility are given by the isospin density and its susceptibility are given by and finally the electric charge density and its susceptibility are To proceed we denote the fermion matrix on the lattice with M and introduce the quark chemical potential µ f in the usual way, i.e. as a constant imaginary vector potential in the temporal direction [56]. We then encounter the following derivatives [57] R (1) where all expectation values are evaluated at vanishing chemical potentials, µ f = 0. It follows from symmetry that n f = R (1) f = 0. The diagonal and off-diagional susceptibilities are then written as   where we used the fact that the fermion matrix is the direct product of the fermion matrices for each flavour. Denoting the degenerate light quarks with = u = d, we find finally that We note that for two degenerate light flavours the isospin susceptibility χ I does not depend on the disconnected term R (3) , while this term contributes more strongly to the baryon susceptibility χ B than to the charge susceptibility χ Q . Note that for three degenerate flavours, ss . The disconnected term is numerically the most expensive quantity to be computed and it dominates the uncertainty of the final results.
We have determined the susceptibilities numerically on our N s = 24 ensembles, see Table 1. The traces in Eq. (4.6) are estimated stochastically, using N v = 9 noise vectors for the connected terms R (2,4) . For the disconnected term R (3) , we use N v = 200 noise vectors in the N τ = 40 case and N v = 100 at the other temperatures. More technical details can be found in Ref. [57].
χ/χ SB Figure 4. Isospin, charge and baryon number susceptibilities (left) and quark number susceptibility for light and strange quarks (right), normalised with the corresponding quantity for free lattice fermions, denoted as χ SB .
In Fig. 3 we present the results for R (2,3,4) , for both light and strange quarks. We observe that at high temperature, the dominant contribution comes from R (2) and that R (3) is compatible with zero within errors, for both the diagonal and off-diagonal components. In this context we note that in hard thermal loop (HTL) perturbation theory [58] the offdiagonal susceptibility is nonzero, showing a clear correlation between different flavours. Also recent lattice calculations [53] have shown a clear dip for the off-diagonal term in the crossover region. Our result might be due to the relatively heavy sea quark masses.
The various susceptibilities χ I , χ Q and χ B are presented in Fig. 4 (left), where all observables are normalised with the corresponding quantities for free lattice fermions, with the same lattice geometry [57]. In the free case the bare fermion anisotropy is set equal to the renormalised value for our ensembles, while the bare quark mass is set to zero. We have evaluated the effect of the uncertainty in the determination of ξ, see Ref. [59], on our final results and found it to be a systematic effect of the order of 5%. In Fig. 4 we clearly see a steep increase above 150 MeV, and for T 250 MeV the value of the susceptibilities is around 85% of the Stefan-Boltzmann value, i.e. the free case. The result for χ Q will be used in Sec. 6 to determine the diffusion coefficient. The baryon number susceptibility behaves qualitatively in a similar way to the other two, but has larger errors due to the way the various terms combine, see Eq. (4.8). Our results appear consistent with the findings of other lattice groups [53,54] and also with resummed perturbation theory [60], in particular concerning the deviation from unity at the highest temperatures.
Finally, in Fig. 4 (right) we show separately the quark number susceptibilities for light and strange quarks, again normalised with the corresponding quantity for free lattice fermions. We see some indication for "flavour separation" during the QCD crossover transition, as discussed in Ref. [61] and reported in Ref. [62], where a continuum extrapolated lattice QCD calculation was performed (see, however, Ref. [63]).

Conserved vector currents and conductivity
We consider the electromagnetic current for three flavours, where j f µ are the vector currents for each flavour and eq f are the corresponding electric charges. The Euclidean current-current correlator G em µν (τ ) is then defined, at zero spatial momentum, as This correlator admits a spectral representation of the form with the kernel where ρ em µν (ω) is the spectral function. Application of linear response theory [13] yields the Kubo formula for the electrical conductivity σ, where the spectral function ρ em (ω) has to be obtained from the Euclidean correlator G em (τ ) = i G em ii (τ ) by inverting Eq. (5.3). It will be useful to normalise the electromagnetic observables by the sum of the square of the individual quark charges, which equals 2e 2 /3 for three flavours. We then define and consider G(τ ) and ρ(ω) from now on. Where the light/strange quark contributions are shown separately, the corresponding correlators and spectral functions are normalised with the electromagnetic prefactor for two light quarks/one strange quark respectively.

Correlators
We use the exactly conserved vector current on the lattice as an interpolator for j em µ , since it is protected from renormalisation under quantum corrections. It is defined as where c 4 = 1/2, c i = 1/(2γ f ) and U µ (x) are the gauge links. The two connected diagrams that contribute to the correlator (5.2) are where S(x, y) = ψ(x)ψ(y) is the fermion propagator, Γ ± µ = 1 ± γ µ , Γ ± µ = 1 ± γ µ , and γ µ = γ 4 γ µ γ 4 . In the following we neglect the disconnected pieces, which is expected to have a small effect, since their contribution is identically zero in the (degenerate) N f = 3 case (since f q f = 0). We note that the same choice has been made in all previous studies [18][19][20][21][22]. Finally, as we have shown in Sec. 4, the contribution from disconnected diagrams to the charge susceptibility is negligible.
We have computed the conserved vector current correlator G(τ ) = i G ii (τ ) for the ensembles presented in Table 2. In Fig. 5  To study the effect of increasing the temperature, we show in Fig. 6 the vector correlators normalised by the free (noninteracting) correlators on the lattice, again for both the light and strange quarks. We observe a clear difference between the low temperature phase, where this ratio decreases with increasing τ , and the high temperature region where the ratio is relatively constant and close to unity, demonstrating that the quarks are quasi-free. The effect of the heavier strange quark mass is clearly visible, both below and above T c . For the light quarks, we observe that all four correlators above T c follow the same pattern and exceed the free correlator by about 7%. This may partly be due to the choice of bare parameters in the free lattice calculation, where we choose the bare anisotropy equal to the renormalised one and the bare quark massm 0 = 0. On the other hand, a similar enhancement above the free correlator has been observed analytically in a next-to-leading order perturbative calculation [64]. For the strange quark, the ratio at the highest temperatures is consistent with 1.
At four temperatures, three above and one below T c , we have access to different spatial volumes, namely L s ∼ 2.9 respectively 3.9 fm, or N s = 24 and 32. To study the finite-size effects, we show in Fig. 7   that finite-size effects are at the percent level or less and decrease at higher temperature.

Spectral functions and conductivity
To obtain the spectral functions and conductivity from the correlators, the spectral representation (5.3) has to be inverted. For this we follow the same procedure as in our previous work [19,23], namely the Maximum Entropy Method (MEM) [65]. Other possibilities include the use of a physically motivated Ansatz for the spectral function, with a number of parameters to be determined [20][21][22], as well as alternative inversion methods [66][67][68].
Since the implementation of MEM has been presented in previous works [19,23,65], we summarise here only the main ingredients. At large ω and nonzero τ , the kernel K(τ, ω) is exponentially suppressed, hence one may impose an upper limit ω < ω max . The finite interval 0 ≤ ω < ω max is then discretised using N ω points. Typical values are a τ ω max = 3 and N ω = 1000. Eq. (5.3) has the form of a generalised Laplace transform, the inverse of which is known to be an ill-posed problem. In MEM one extracts the most probable spectral function ρ(ω), given some prior knowledge H and the data D. This is expressed where L = 1 2 χ 2 is the standard likelihood function and S is the Shannon-Jaynes entropy, Here m(ω) is the default model which implements the prior information on ρ(ω) in the absence of the data D. The result for ρ(ω) is then obtained by extremising P [ρ|DH]. To do this, we use a modification [19] of Bryan's algorithm [69] which cures the 1/ω instability of the kernel K(τ, ω) at small ω. The default model we use is [19] m(ω) = m 0 (b + ω)ω, (5.12) where m 0 is an overall normalisation set by fitting the correlator to the trial function obtained by using m(ω) in the convolution integral (5.3). This default model is chosen because it matches the perturbative ω-dependence at large ω in the continuum theory (on the lattice this behaviour is modified due to the finite Brillouin zone [70,71]) and allows a nonzero value of ρ(ω)/ω as ω → 0 and hence a nonzero conductivity σ according to the Kubo relation (5.5). As always, it is essential to check that the spectral functions obtained are independent of the choices made in the MEM procedure, including the choice of default model and its parameters. This is discussed in some detail in the next subsection, after presenting the results. The spectral functions obtained with MEM are shown in Fig. 8, normalised as ρ(ω)/ω 2 . The spectral functions are shown for light (left) and strange (right) quarks at three representative temperatures spanning the entire range. We always use the largest volume, N s = 32, when available. The vertical dashed lines correspond to an estimate of the mass of the ground state in the corresponding vector channel at zero temperature [39]. The MEM analysis indeed indicates a peak at this value below T c , which becomes less pronounced and disappears as the temperature increases. The divergence at small ω at the higher temperatures is due to the transport peak. This is emphasised in Fig. 9 where ρ(ω)/ωT is shown. According to the Kubo relation (5.5), the intercepts are proportional to σ/T . We observe a conductivity which is clearly nonzero above T c and which depends on the quark mass.
The final results for the conductivity are shown in Fig. 10 as a function of the temperature. We present the result as C −1 em σ/T where C em was defined in Eq. (5.6). The results are shown for the light and strange quarks separately and for all three quarks combined. Note that we always first construct the electromagnetic current operator with the correct T /T c Figure 10. Temperature dependence of C −1 em σ/T for light and strange quarks separately, slightly shifted for clarity (above) and combined (below). The vertical size of the rectangles reflects the systematic uncertainty due to changes in the default model, while the whiskers depict the statistical jackknife error on top of this.
weighting of the quark charges and then apply MEM to the resulting correlators. The systematic uncertainty due to the choice of the parameter b in the default model (5.12) is represented by the vertical size of the filled rectangles. This is discussed further below. The statistical uncertainty due to the finite number of configurations is represented by the upper and lower whiskers emanating from the rectangles, so that the total spread of values (statistical and systematic) is given by the size of the error bars.
We observe that the contributions from the light and the strange quarks are comparable, except in the crossover region, where the strange quark contribution is suppressed. Note, however, that in the total conductivity the strange quark contribution is in any case suppressed with respect to the light quarks, due to the different electromagnetic prefactors: C em = 5e 2 /9 versus C s em = e 2 /9.

Systematics
In order to have confidence in the results, it is necessary to study the systematic uncertainties in the MEM analysis. As in our previous work, the sensitivity to ω max is modest, provided that 3 a τ ω max 5. Here we show results from tests varying the b parameter in the default model, excluding intermediate time points, and varying the choice of time range included in the analysis.
In Fig. 11 we the systematic error coming from the default model in our conductivity determination, see Fig. 10. On our anisotropic lattices, the temporal lattice spacing is smaller than the spatial one, with ξ = a s /a τ = 3.5. Hence we have more time slices available for the MEM analysis than on an isotropic lattice with the same spatial lattice spacing. One may question how this improves the results, if at all. We test this by including in the MEM analysis either all time slices (always discarding the first four), or one in two, or one in three. With our specific value of ξ, the latter is roughly equivalent to the isotropic case. The results are shown in Fig. 12. We observe that at the lower temperatures the results are manifestly stable and consistent and hence an isotropic lattice would suffice. On the other hand, at the higher temperatures the benefit of the anisotropy is clearly visible: while not affecting the central value substantially, it greatly reduces the systematic uncertainty in the reconstruction. These results indicate the robustness of the results and the necessity of using anisotropic lattices.
Finally, we assess the uncertainty arising from the choice of time window used in the MEM analysis. Since we work at a fixed lattice spacing, increasing the temperature implies having fewer time slices available. Hence, one possibility could be that the observed temperature dependence of the conductivity is simply an artefact due to the different number of time slices available and hence not physical. We test this by using MEM with restricted time windows at various temperatures, see Fig. 13. In each row we consider three

Diffusion coefficient
We are now in the position to combine the results for the conductivity σ and the charge susceptibility χ Q to obtain the charge diffusion coefficient D = σ/χ Q . This ratio is independent of the electromagnetic factor C em . We present the result for the dimensionless combination 2πT D in Fig. 14. The vertical size of the rectangles represents the systematic uncertainty coming from the determination of the conductivity, while the whiskers indicate the statistical uncertainty in both σ and χ Q . The first observation we make is that the diffusion coefficient is of the order of 1/(2πT ). In order to judge whether this is a sensible result, we note that in strongly coupled theories, in particular those which can be treated with holography, this is exactly the magnitude that is expected. For instance, the diffusion coefficient for R-charge in N = 4 Yang-Mills theory at nonzero temperature equals 1/(2πT ) [32][33][34]. On the other hand, in weakly coupled theories the diffusion coefficient, being proportional to the mean free path, is large and diverges as the interactions are turned off. Hence our results are consistent with strongly coupled real-time dynamics in the quark-gluon plasma in this temperature range and in the hydrodynamic regime.
The second observation is a dip in the diffusion coefficient in the transition region. This can be understood as follows. We first note that at high temperature DT is expected to rise, since the conductivity is expected to increase, due to the diverging mean free path at weak coupling, while the susceptibility is of the order of the Stefan-Boltzmann value. Their ratio will hence grow large. On the low-temperature side we note that the susceptibility drops rapidly in the confined phase. On the other hand, we expect the conductivity to be nonzero, since it can be assumed that a pion gas at low temperature is a conductor rather than an insulator. Hence the ratio will again lead to a rise of DT . This then naturally leads to a minimum around T c , as in the case of the ratio of the shear viscosity to entropy density [3,72]. As a side remark we note that a successful numerical evaluation at very low temperature along the lines followed here will be unlikely, since D involves the ratio of two suppressed quantities in the confined phase.
We note here that a plot similar to Fig. 14 was constructed in Ref. [4], by combining the conductivity results for the two light flavours from Ref. [23] with the (continuumextrapolated) susceptibility results from Ref. [53]. The conductivity and (quark number) susceptibility have also been computed in Ref. [20] for quenched QCD and in Ref. [22] for QCD with N f = 2 flavours, but the resulting diffusion coefficient was not given. Note that the latter also contains a comparison with N = 4 Yang-Mills theory. T /T c Figure 14. Diffusion coefficient D multiplied by 2πT as a function of the temperature T , using D = σ/χ Q . The vertical size of the rectangles represents the systematic uncertainty due to the uncertainty in the estimate of the conductivity (see Fig. 10), while the whiskers indicate the statistical jackknife error from both σ and χ Q .
Finally, we remark that an attempt to determine the charm diffusion coefficient can be found in Ref. [73] using quenched lattice simulations on large and fine isotropic lattices, with the finding that D ∼ 1/(πT ) in the deconfined phase. For very heavy quarks various diffusion coefficients are being determined using heavy-quark effective theory, see e.g. Refs. [74][75][76] and references therein.

Conclusion
The main result in this paper is the determination of the electrical conductivity and charge diffusion coefficient at nonzero temperature in QCD with N f = 2 + 1 quark flavours, using anisotropic lattice QCD simulations.
Our results for the conductivity σ confirm our previous findings where only the u and d quark contributions were taken into account: σ/T increases by a factor of 5-6 in our temperature range, which spans the chiral and deconfinement transition. We note that the results for the conductivity at the lowest temperature should be treated with caution, since a possible narrow transport peak resulting from hadronic interactions would not be detectable with our methods. We find that the diffusion coefficient is of the order of 1/(2πT ) and has a dip around the transition temperature between the confined and the deconfined phase. This is consistent with a strongly-coupled quark-gluon plasma in the hydrodynamic regime.
In order to reach this result, we have used the Maximum Entropy Method to construct spectral functions from the numerically determined Euclidean correlators of the conserved vector current. The conductivity then follows from the linear behaviour of the spectral functions at small energies. Independently we have determined various second-order susceptibilities and found agreement with previous results. The diffusion coefficient is given by the ratio of the electrical conductivity to the charge susceptibility.
As an outlook, we note that there are various things that can be improved. Besides MEM, it might be useful to apply other recently developed inversion methods [66][67][68]. It should be stated that our results are robust against variation of several systematic input variables, most notably those related to the Euclidean time interval and number of Euclidean time points included. Here we found that the anisotropy is essential at the highest temperatures. Concerning our ensembles, we note that the spatial lattice is relatively coarse and that the light quarks are heavier than in nature, which affects the transition temperature. Hence it is worthwhile to repeat the analysis with lighter quarks on finer lattices. An increase of the anisotropy on the other hand will allow for even better control on systematics of the inversion. We hope to address some of these issues in the future.