# Electrical conductivity and charge diffusion in thermal QCD from the lattice

###### Abstract

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 and has a dip around the crossover temperature. Our results are obtained with lattice simulations containing dynamical flavours on anisotropic lattices. The Maximum Entropy Method is used to construct spectral functions from correlators of the conserved vector current.

###### Keywords:

## 1 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 Heinz:2013th (); Gale:2013da (). 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 is close to the value obtained in models with holographic duals Kovtun:2004de ().

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 Ling:2013ksb () and the evolution of strong electromagnetic fields produced in noncentral collisions Tuchin:2013ie (); McLerran:2013hla (); Gursoy:2014aka (); Zakharov:2014dia (); Tuchin:2014iua (); Satow:2014lia (). It has also been suggested that experimental information on conductivity can be extracted from flow parameters in heavy-ion collisions Hirono:2012rt ().

Using linear response theory, transport coefficients can be related to current-current spectral functions in thermal equilibrium, via the celebrated Kubo relations kubo (); Kadanoff1963419 (). 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 Aarts:2002cc (). In recent years several results have been obtained along these lines. Refs. Meyer:2007ic (); Meyer:2007dy () contain the first (theoretically controlled) results for the shear and bulk viscosity, see also the review Meyer:2011gj (). 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 Gupta:2003zh (); Aarts:2007wj () in quenched QCD. Since then the study of the systematic uncertainties and the extension to Wilson-type quarks have taken a central role, in quenched QCD Ding:2010ga (); Kaczmarek:2013dya () and in dynamical QCD with and flavours Brandt:2012jc (); Amato:2013naa (). All studies Aarts:2007wj (); Ding:2010ga (); Kaczmarek:2013dya (); Brandt:2012jc (); Amato:2013naa () are in qualitative agreement: around , where is the crossover temperature, the conductivity is of the order of , where is an electromagnetic prefactor depending on the quark charges (see Sec. 5 below). Recent non-lattice studies include Refs. Cassing:2013iz (); Finazzo:2013efa (); Steinert:2013fza (); Yin:2013kya (); Lee:2014pwa (); Puglisi:2014pda (); Puglisi:2014sha (); Greif:2014oia ().

The most detailed lattice study so far can be found in Ref. Amato:2013naa (), where the temperature dependence of across the deconfinement transition was studied for the first time, over a range of temperatures corresponding to . In that analysis we used a QGP with 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 , which allows us to compute for the first time the charge diffusion coefficient 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 in Yang-Mills theory at nonzero temperature Policastro:2002se (); Teaney:2006nc (); Son:2007vk (), 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.

## 2 lattice details

In this work we follow the Hadron Spectrum Collaboration Edwards:2008ja () and use a Symanzik-improved anisotropic gauge action with tree-level mean-field coefficients and a mean-field–improved Wilson-clover fermion action with stout-smeared links Morningstar:2003gk (). The anisotropy, with a reduced temporal lattice spacing , 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, , common to the gauge and fermionic degrees of freedom. The ensembles employed here are part of our “2nd generation” data set Allton:2014uia () and were previously used for a determination of the conductivity (from two light flavours only) Amato:2013naa () and the bottomonium spectrum at nonzero temperature Aarts:2014cda ().

The gauge action takes the form

(2.1) |

where and are the plaquette and rectangular Wilson loops, are the spatial (temporal) mean links, is the bare gauge anisotropy and, as usual, with colours.

The fermion action (for a single flavour) reads

(2.2) |

where is the bare fermion mass, the bare fermion anisotropy, covariant finite differences, , and the clover coefficients

(2.3) |

The spatial gauge links in the fermion action have been stout smeared Morningstar:2003gk () with smearing weight and iterations, and is the mean value of the spatial stout-smeared links.

gauge coupling | |||

bare gauge anisotropy | = fm | ||

bare fermion anisotropy | GeV | ||

spatial clover coefficient | = 3.5 | ||

temporal clover coefficient | MeV | ||

bare light quark mass | |||

bare strange quark mass |

The choice of bare parameters is given in Table 1 and follows the tuning by the Hadron Spectrum Collaboration Edwards:2008ja (). The resulting renormalised anisotropy is . The two degenerate light quarks yield MeV (2.8 times larger than the physical pion), corresponding to Lin:2008pr (), while the third flavour is tuned to the strange quark mass Edwards:2008ja ().

[MeV] | ||||||
---|---|---|---|---|---|---|

24 | 32 | 16 | 352 | 1.90 | 1059 | 4 |

24 | 20 | 281 | 1.52 | 1001 | 4 | |

24 | 32 | 24 | 235 | 1.27 | 500 | 4 |

24 | 32 | 28 | 201 | 1.09 | 502 | 4 |

24 | 32 | 32 | 176 | 0.95 | 501 | 4 |

24 | 36 | 156 | 0.84 | 501 | 4 | |

24 | 40 | 141 | 0.76 | 523 | 4 | |

32 | 48 | 117 | 0.63 | 601 | 4 | |

24 | 128 | 44 | 0.24 | 401 | 1 |

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 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 respectively 3.9 fm ( and ). Access to the zero-temperature configurations () 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 Edwards:2004sx () with Bagel routines Boyle:2009vp (). For further details about the algorithm, we refer to Refs. Edwards:2008ja (); Lin:2008pr (). After 1000 thermalisation trajectories (2000 for the ensemble), configurations were sampled every 10 RHMC trajectories, except for the 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.

## 3 Deconfinement and chiral transition

After a determination of the lattice spacing, the temperature can be specified in MeV very precisely using the standard relation . However, it is desirable to express the temperature in units of , 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. Aoki:2006we ()). 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. Borsanyi:2012xf ().

The Polyakov loop expectation value is related to the free energy of a static quark via

(3.1) |

However, is only defined up to an additive renormalisation constant , which depends on the gauge coupling and other bare parameters but not on the temperature. Expressing the renormalised free energy as allows us write the renormalised Polyakov loop as

(3.2) |

which defines the multiplicative renormalisation constant . Following Ref. Borsanyi:2012xf (), we impose a renormalisation condition at a reference temperature , by requiring that

(3.3) |

which fixes .

Figure 1 shows the Polyakov loop with three different renormalisation schemes corresponding to different choices of 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 or , where the error reflects the systematic error coming from the spread of the three renormalisation schemes. This corresponds to a deconfinement critical temperature of 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 are commonly used to define the crossover transition temperature. We have determined the chiral susceptibility due to the two light flavours, using Aoki:2006br (); Borsanyi:2010bp (); Bazavov:2011nk ()

(3.4) | |||||

where is the partition function, the fermion matrix, the spatial volume and the degenerate light quark mass. Moreover, we introduce here the connected and the disconnected contributions to the susceptibility. The traces in Eq. (3.4) are determined using 16 noise vectors for the disconnected contribution and 4 for the connected one. Because we change the temperature by changing the value of rather than the lattice spacing, the (additive and multiplicative) renormalisation of is the same for all temperatures. A peak in the susceptibility occurs therefore at the same temperature for the renormalised and unrenormalised . 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 MeV, which is somewhat lower than the value obtained from the Polyakov loop. In the following we will use the value for determined from the Polyakov loop as this has been obtained with a better precision.

## 4 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 Asakawa:2000wh () by studying event-by-event fluctuations of charged particle ratios Bleicher:2000ek (). 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 Gottlieb:1988cq (); Allton:2002zi (); Gavai:2005sd (); Bernard:2007nm (); Borsanyi:2011sw (); Bazavov:2013uja (). 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. Borsanyi:2012uq (). The charge diffusion coefficient and the electrical conductivity are related via the well-known relation , where is the charge susceptibility Kadanoff1963419 (). In this section we determine 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

(4.1) |

where is the partition function, the spatial volume, and the quark chemical potentials for flavours . Baryon (), isospin () and electrical charge () chemical potentials are related to the quark chemical potentials as

(4.2) |

In general we denote the electrical charge of the quark as , with the elementary charge and or its fractional charge.

All desired quantities can now be expressed in terms of and . The baryon number density and baryon number susceptibility are given by

(4.3) |

the isospin density and its susceptibility are given by

(4.4) |

and finally the electric charge density and its susceptibility are

(4.5) |

To proceed we denote the fermion matrix on the lattice with and introduce the quark chemical potential in the usual way, i.e. as a constant imaginary vector potential in the temporal direction Hasenfratz:1983ba (). We then encounter the following derivatives Giudice:2013fza ()

(4.6) |

where all expectation values are evaluated at vanishing chemical potentials, . It follows from symmetry that . The diagonal and off-diagional susceptibilities are then written as

(4.7) |

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 , we find finally that

(4.8) |

We note that for two degenerate light flavours the isospin susceptibility does not depend on the disconnected term , while this term contributes more strongly to the baryon susceptibility than to the charge susceptibility . Note that for three degenerate flavours, is also independent of , since . 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 ensembles, see Table 1. The traces in Eq. (4.6) are estimated stochastically, using noise vectors for the connected terms . For the disconnected term , we use noise vectors in the case and at the other temperatures. More technical details can be found in Ref. Giudice:2013fza ().

In Fig. 3 we present the results for , for both light and strange quarks. We observe that at high temperature, the dominant contribution comes from and that 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 Blaizot:2001vr () the off-diagonal susceptibility is nonzero, showing a clear correlation between different flavours. Also recent lattice calculations Borsanyi:2011sw () 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 , and are presented in Fig. 4 (left), where all observables are normalised with the corresponding quantities for free lattice fermions, with the same lattice geometry Giudice:2013fza (). 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. Dudek:2012xn (), on our final results and found it to be a systematic effect of the order of . In Fig. 4 we clearly see a steep increase above MeV, and for MeV the value of the susceptibilities is around of the Stefan-Boltzmann value, i.e. the free case. The result for 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 Borsanyi:2011sw (); Bazavov:2013uja () and also with resummed perturbation theory Andersen:2012wr (), 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. Ratti:2011au () and reported in Ref. Bellwied:2013cta (), where a continuum extrapolated lattice QCD calculation was performed (see, however, Ref. Bazavov:2013dta ()).

## 5 Conserved vector currents and conductivity

We consider the electromagnetic current for three flavours,

(5.1) |

where are the vector currents for each flavour and are the corresponding electric charges. The Euclidean current-current correlator is then defined, at zero spatial momentum, as

(5.2) |

This correlator admits a spectral representation of the form

(5.3) |

with the kernel

(5.4) |

where is the spectral function.

Application of linear response theory Kadanoff1963419 () yields the Kubo formula for the electrical conductivity ,

(5.5) |

where the spectral function has to be obtained from the Euclidean correlator 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,

(5.6) |

which equals for three flavours. We then define

(5.7) |

and consider 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.

### 5.1 Correlators

We use the exactly conserved vector current on the lattice as an interpolator for , since it is protected from renormalisation under quantum corrections. It is defined as

(5.8) |

where , and are the gauge links. The two connected diagrams that contribute to the correlator (5.2) are

(5.9) |

where is the fermion propagator, , , and . 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) case (since ). We note that the same choice has been made in all previous studies Gupta:2003zh (); Aarts:2007wj (); Ding:2010ga (); Brandt:2012jc (); Kaczmarek:2013dya (). 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 for the ensembles presented in Table 2. In Fig. 5 the results are shown for the light quarks (left) and strange quark (right) at all the temperatures, as a function of Euclidean time. Note that we present the correlators at the different temperatures as a function of , which has the effect of separating them, even when they have identical decay. The correlators are symmetric about .

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 . For the light quarks, we observe that all four correlators above 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 mass . On the other hand, a similar enhancement above the free correlator has been observed analytically in a next-to-leading order perturbative calculation Burnier:2014cna (). For the strange quark, the ratio at the highest temperatures is consistent with 1.

At four temperatures, three above and one below , we have access to different spatial volumes, namely 2.9 respectively 3.9 fm, or and 32. To study the finite-size effects, we show in Fig. 7 the ratio of the to the correlators. We observe that finite-size effects are at the percent level or less and decrease at higher temperature.

### 5.2 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 Aarts:2007wj (); Amato:2013naa (), namely the Maximum Entropy Method (MEM) mem (). Other possibilities include the use of a physically motivated Ansatz for the spectral function, with a number of parameters to be determined Brandt:2012jc (); Ding:2010ga (); Kaczmarek:2013dya (), as well as alternative inversion methods Burnier:2011jq (); Burnier:2012ts (); Burnier:2013nla ().

Since the implementation of MEM has been presented in previous works mem (); Aarts:2007wj (); Amato:2013naa (), we summarise here only the main ingredients. At large and nonzero , the kernel is exponentially suppressed, hence one may impose an upper limit . The finite interval is then discretised using points. Typical values are and . 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 and the data . This is expressed as a conditional probability, via Bayes theorem,

(5.10) |

where is the standard likelihood function and is the Shannon-Jaynes entropy,

(5.11) |

Here is the default model which implements the prior information on in the absence of the data . The result for is then obtained by extremising . To do this, we use a modification Aarts:2007wj () of Bryan’s algorithm bryan () which cures the instability of the kernel at small . The default model we use is Aarts:2007wj ()

(5.12) |

where is an overall normalisation set by fitting the correlator to the trial function obtained by using 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 Karsch:2003wy (); Aarts:2005hg ()) and allows a nonzero value of as 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 . 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, , 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 Lin:2008pr (). The MEM analysis indeed indicates a peak at this value below , 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 is shown. According to the Kubo relation (5.5), the intercepts are proportional to . We observe a conductivity which is clearly nonzero above 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 where 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 weighting of the quark charges and then apply MEM to the resulting correlators. The systematic uncertainty due to the choice of the parameter 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: versus .

### 5.3 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 is modest, provided that . Here we show results from tests varying the parameter in the default model, excluding intermediate time points, and varying the choice of time range included in the analysis.

In Fig. 11 we show the dependence of the conductivity on the default model parameter . The results are stable provided . We therefore use the range to define 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 . 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 temperatures. We always perform the MEM analysis starting at . We then include either all time slices available (filled and open symbols in the left plots), or constrain the number of time slices by the highest temperature in each row (filled symbols only). In this way we can study the effect of adding more time points as the temperature is decreased. The resulting conductivities are shown on the right. Here the filled symbols are obtained using the restricted MEM analysis, while the open symbols indicate the results with all available time slices used. The results are remarkably stable, with agreement between open and filled symbols at each temperature, including the highest one shown. We hence conclude that the observed dependence of the conductivity on the temperature is a thermal effect, rather than a bias introduced in our method.

## 6 Diffusion coefficient

We are now in the position to combine the results for the conductivity and the charge susceptibility to obtain the charge diffusion coefficient . This ratio is independent of the electromagnetic factor . We present the result for the dimensionless combination 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 .

The first observation we make is that the diffusion coefficient is of the order of . 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 -charge in Yang-Mills theory at nonzero temperature equals ) Policastro:2002se (); Teaney:2006nc (); Son:2007vk (). 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 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 . This then naturally leads to a minimum around , as in the case of the ratio of the shear viscosity to entropy density Kovtun:2004de (); Csernai:2006zz (). As a side remark we note that a successful numerical evaluation at very low temperature along the lines followed here will be unlikely, since 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. Ling:2013ksb (), by combining the conductivity results for the two light flavours from Ref. Amato:2013naa () with the (continuum-extrapolated) susceptibility results from Ref. Borsanyi:2011sw (). The conductivity and (quark number) susceptibility have also been computed in Ref. Ding:2010ga () for quenched QCD and in Ref. Brandt:2012jc () for QCD with flavours, but the resulting diffusion coefficient was not given. Note that the latter also contains a comparison with Yang-Mills theory.

Finally, we remark that an attempt to determine the charm diffusion coefficient can be found in Ref. Ding:2012sp () using quenched lattice simulations on large and fine isotropic lattices, with the finding that in the deconfined phase. For very heavy quarks various diffusion coefficients are being determined using heavy-quark effective theory, see e.g. Refs. CaronHuot:2009uh (); Banerjee:2011ra (); Francis:2013cva () and references therein.

## 7 Conclusion

The main result in this paper is the determination of the electrical conductivity and charge diffusion coefficient at nonzero temperature in QCD with quark flavours, using anisotropic lattice QCD simulations.

Our results for the conductivity confirm our previous findings where only the and quark contributions were taken into account: 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 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 Burnier:2011jq (); Burnier:2012ts (); Burnier:2013nla (). 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.

## Acknowledgements

We thank Benjamin Jäger and Prem Kumar for discussion. We are grateful to the Hadron Spectrum Collaboration for providing the zero-temperature ensemble. For computational resources we thank the STFC funded DiRAC Facility, HPC Wales, PRACE via grants 2011040469 and 2012061129, and ICHEC. GA is supported by STFC, the Royal Society, the Wolfson Foundation and the Leverhulme Trust. CRA is supported by STFC and the Leverhulme Trust. He thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work. AA was supported by a Swansea University postgraduate scholarship and acknowledges European Union Grant Agreement number 238353 (ITN STRONGnet) and Academy of Finland grant 1267286 for support during the completion of this work. SH is supported by STFC.

## References

- (1) U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. 63 (2013) 123 [arXiv:1301.2826 [nucl-th]].
- (2) C. Gale, S. Jeon and B. Schenke, Int. J. Mod. Phys. A 28 (2013) 1340011 [arXiv:1301.5893 [nucl-th]].
- (3) P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94 (2005) 111601 [hep-th/0405231].
- (4) B. Ling, T. Springer and M. Stephanov, Phys. Rev. C 89 (2014) 064901 [arXiv:1310.6036 [nucl-th]].
- (5) K. Tuchin, Adv. High Energy Phys. 2013 (2013) 490495 [arXiv:1301.0099].
- (6) L. McLerran and V. Skokov, Nucl. Phys. A 929 (2014) 184 [arXiv:1305.0774 [hep-ph]].
- (7) U. Gursoy, D. Kharzeev and K. Rajagopal, Phys. Rev. C 89 (2014) 054905 [arXiv:1401.3805 [hep-ph]].
- (8) B. G. Zakharov, Phys. Lett. B 737 (2014) 262 [arXiv:1404.5047 [hep-ph]].
- (9) K. Tuchin, arXiv:1411.1363 [hep-ph].
- (10) D. Satow, Phys. Rev. D 90 (2014) 034018 [arXiv:1406.7032 [hep-ph]].
- (11) Y. Hirono, M. Hongo and T. Hirano, Phys. Rev. C 90 (2014) 2, 021903 [arXiv:1211.1114 [nucl-th]].
- (12) R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
- (13) L. P. Kadanoff and P. C. Martin, Annals of Physics 24, 419 (1963).
- (14) G. Aarts and J. M. Martínez Resco, JHEP 0204 (2002) 053 [hep-ph/0203177].
- (15) H. B. Meyer, Phys. Rev. D 76 (2007) 101701 [arXiv:0704.1801 [hep-lat]].
- (16) H. B. Meyer, Phys. Rev. Lett. 100 (2008) 162001 [arXiv:0710.3717 [hep-lat]].
- (17) H. B. Meyer, Eur. Phys. J. A 47 (2011) 86 [arXiv:1104.3708 [hep-lat]].
- (18) S. Gupta, Phys. Lett. B 597 (2004) 57 [hep-lat/0301006].
- (19) G. Aarts, C. Allton, J. Foley, S. Hands and S. Kim, Phys. Rev. Lett. 99 (2007) 022002 [hep-lat/0703008].
- (20) H.-T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann and W. Soeldner, Phys. Rev. D 83 (2011) 034504 [arXiv:1012.4963 [hep-lat]].
- (21) O. Kaczmarek and M. Müller, PoS LATTICE 2013 (2013) 175 [arXiv:1312.5609 [hep-lat]].
- (22) B. B. Brandt, A. Francis, H. B. Meyer and H. Wittig, JHEP 1303 (2013) 100 [arXiv:1212.4200 [hep-lat]].
- (23) A. Amato, G. Aarts, C. Allton, P. Giudice, S. Hands and J. I. Skullerud, Phys. Rev. Lett. 111 (2013) 172001 [arXiv:1307.6763 [hep-lat]].
- (24) W. Cassing, O. Linnyk, T. Steinert and V. Ozvenchuk, Phys. Rev. Lett. 110 (2013) 18, 182301 [arXiv:1302.0906 [hep-ph]].
- (25) S. I. Finazzo and J. Noronha, Phys. Rev. D 89 (2014) 106008 [arXiv:1311.6675 [hep-th]].
- (26) T. Steinert and W. Cassing, Phys. Rev. C 89 (2014) 035203 [arXiv:1312.3189 [hep-ph]].
- (27) Y. Yin, Phys. Rev. C 90 (2014) 044903 [arXiv:1312.4434 [nucl-th]].
- (28) C. H. Lee and I. Zahed, Phys. Rev. C 90 (2014) 025204 [arXiv:1403.1632 [hep-ph]].
- (29) A. Puglisi, S. Plumari and V. Greco, arXiv:1407.2559 [hep-ph].
- (30) A. Puglisi, S. Plumari and V. Greco, Phys. Rev. D 90 (2014) 114009 [arXiv:1408.7043 [hep-ph]].
- (31) M. Greif, I. Bouras, C. Greiner and Z. Xu, Phys. Rev. D 90 (2014) 9, 094014 [arXiv:1408.7049 [nucl-th]].
- (32) G. Policastro, D. T. Son and A. O. Starinets, JHEP 0209 (2002) 043 [hep-th/0205052].
- (33) D. Teaney, Phys. Rev. D 74 (2006) 045025 [hep-ph/0602044].
- (34) D. T. Son and A. O. Starinets, Ann. Rev. Nucl. Part. Sci. 57 (2007) 95 [arXiv:0704.0240 [hep-th]].
- (35) R. G. Edwards, B. Joo and H. W. Lin, Phys. Rev. D 78 (2008) 054501 [arXiv:0803.3960 [hep-lat]].
- (36) C. Morningstar and M. J. Peardon, Phys. Rev. D 69 (2004) 054501 [hep-lat/0311018].
- (37) C. Allton, G. Aarts, A. Amato, W. Evans, P. Giudice, T. Harris, S. Hands and A. Kelly et al., PoS LATTICE 2013 (2014) 151 [arXiv:1401.2116 [hep-lat]].
- (38) G. Aarts, C. Allton, T. Harris, S. Kim, M. P. Lombardo, S. M. Ryan and J. I. Skullerud, JHEP 1407 (2014) 097 [arXiv:1402.6210 [hep-lat]].
- (39) H. W. Lin et al. [Hadron Spectrum Collaboration], Phys. Rev. D 79 (2009) 034502 [arXiv:0810.3588 [hep-lat]].
- (40) R. G. Edwards et al. [SciDAC and LHPC and UKQCD Collaborations], Nucl. Phys. Proc. Suppl. 140 (2005) 832 [hep-lat/0409003].
- (41) P. A. Boyle, Comput. Phys. Commun. 180 (2009) 2739.
- (42) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443 (2006) 675 [hep-lat/0611014].
- (43) S. Borsanyi, Y. Delgado, S. Durr, Z. Fodor, S. D. Katz, S. Krieg, T. Lippert and D. Nogradi et al., Phys. Lett. B 713 (2012) 342 [arXiv:1204.4089 [hep-lat]].
- (44) Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B 643 (2006) 46 [hep-lat/0609068].
- (45) S. Borsanyi et al. [Wuppertal-Budapest Collaboration], JHEP 1009 (2010) 073 [arXiv:1005.3508 [hep-lat]].
- (46) A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta and P. Hegde et al., Phys. Rev. D 85 (2012) 054503 [arXiv:1111.1710 [hep-lat]].
- (47) M. Asakawa, U. W. Heinz and B. Muller, Phys. Rev. Lett. 85 (2000) 2072 [hep-ph/0003169].
- (48) M. Bleicher, S. Jeon and V. Koch, Phys. Rev. C 62 (2000) 061902 [hep-ph/0006201].
- (49) S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken and R. L. Sugar, Phys. Rev. D 38 (1988) 2888.
- (50) C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, C. Schmidt and L. Scorzato, Phys. Rev. D 66 (2002) 074507 [hep-lat/0204010].
- (51) R. V. Gavai and S. Gupta, Phys. Rev. D 72 (2005) 054006 [hep-lat/0507023].
- (52) C. Bernard, C. E. DeTar, L. Levkova, S. Gottlieb, U. M. Heller, J. E. Hetrick, R. Sugar and D. Toussaint, Phys. Rev. D 77 (2008) 014503 [arXiv:0710.1330 [hep-lat]].
- (53) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti and K. Szabo, JHEP 1201 (2012) 138 [arXiv:1112.4416 [hep-lat]].
- (54) A. Bazavov, H.-T. Ding, P. Hegde, F. Karsch, C. Miao, S. Mukherjee, P. Petreczky and C. Schmidt et al., Phys. Rev. D 88 (2013) 9, 094021 [arXiv:1309.2317 [hep-lat]].
- (55) S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, D. Nogradi and K. K. Szabo et al., JHEP 1208 (2012) 126 [arXiv:1205.0440 [hep-lat]].
- (56) P. Hasenfratz and F. Karsch, Phys. Lett. B 125 (1983) 308.
- (57) P. Giudice, G. Aarts, C. Allton, A. Amato, S. Hands and J. I. Skullerud, PoS LATTICE 2013 (2014) 492 [arXiv:1309.6253 [hep-lat]].
- (58) J. P. Blaizot, E. Iancu and A. Rebhan, Phys. Lett. B 523 (2001) 143 [hep-ph/0110369].
- (59) J. J. Dudek, R. G. Edwards and C. E. Thomas, Phys. Rev. D 87 (2013) 3, 034505 [arXiv:1212.0830 [hep-ph]].
- (60) J. O. Andersen, S. Mogliacci, N. Su and A. Vuorinen, Phys. Rev. D 87 (2013) 074003 [arXiv:1210.0912 [hep-ph]].
- (61) C. Ratti, R. Bellwied, M. Cristoforetti and M. Barbaro, Phys. Rev. D 85 (2012) 014004 [arXiv:1109.6243 [hep-ph]].
- (62) R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz and C. Ratti, Phys. Rev. Lett. 111 (2013) 202302 [arXiv:1305.6297 [hep-lat]].
- (63) A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa and S. Mukherjee et al., Phys. Rev. Lett. 111 (2013) 082301 [arXiv:1304.7220 [hep-lat]].
- (64) Y. Burnier, arXiv:1410.1304 [hep-ph].
- (65) M. Asakawa, T. Hatsuda, and Y. Nakahara, Prog. Part. Nucl. Phys. 46 (2001) 459 (2001) [hep-lat/0011040].
- (66) Y. Burnier, M. Laine and L. Mether, Eur. Phys. J. C 71 (2011) 1619 [arXiv:1101.5534 [hep-lat]].
- (67) Y. Burnier and M. Laine, Eur. Phys. J. C 72 (2012) 1902 [arXiv:1201.1994 [hep-lat]].
- (68) Y. Burnier and A. Rothkopf, Phys. Rev. Lett. 111 (2013) 182003 [arXiv:1307.6106 [hep-lat]].
- (69) R.K. Bryan, Eur. Biophys. J. 18 (1990) 165.
- (70) F. Karsch, E. Laermann, P. Petreczky and S. Stickan, Phys. Rev. D 68 (2003) 014504 [hep-lat/0303017].
- (71) G. Aarts and J. M. Martínez Resco, Nucl. Phys. B 726 (2005) 93 [hep-lat/0507004].
- (72) L. P. Csernai, J. I. Kapusta and L. D. McLerran, Phys. Rev. Lett. 97 (2006) 152303 [nucl-th/0604032].
- (73) H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz and W. Soeldner, Phys. Rev. D 86 (2012) 014509 [arXiv:1204.4945 [hep-lat]].
- (74) S. Caron-Huot, M. Laine and G. D. Moore, JHEP 0904 (2009) 053 [arXiv:0901.1195 [hep-lat]].
- (75) D. Banerjee, S. Datta, R. Gavai and P. Majumdar, Phys. Rev. D 85 (2012) 014510 [arXiv:1109.5738 [hep-lat]].
- (76) A. Francis, O. Kaczmarek, M. Laine, M. Müller, T. Neuhaus and H. Ohno, PoS LATTICE 2013 (2013) 453 [arXiv:1311.3759 [hep-lat]].