{'fraction_non_alphanumeric': 0.03143060062109206, 'fraction_numerical': 0.013480204431199431, 'mean_word_length': 4.631715159595326, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Status of Lattice QCD at Finite Temperature

The status of lattice QCD investigations at high temperature is reviewed. After a short introduction into thermal QCD on the lattice we report on the present understanding of the phase diagram and the equation of state, in particular in presence of dynamical quarks. We continue with a discussion of various screening lengths in the plasma phase including results from dimensionally reduced QCD. This is followed by summarizing lattice data on quark number susceptibilities and spectral densities, both of which are of immediate relevance to the interpretation of heavy ion experiments. A major section is devoted to presenting simulations of QCD at small yet phenomenologically important values for the baryon density.


Introduction
Quantumchromodynamics (QCD) is the theory describing the strong interactions, carried by gluons, which confine quarks into hadrons. It is the fundamental theory underlying all nuclear physics, and as such responsible for many salient features of matter as we know it today. A key property of the theory is asymptotic freedom, according to which the coupling strength decreases with the energy transfer of an interaction. One consequence is that at high enough energies of a few GeV the theory is tractable by perturbation theory, which leads to its experimental verification in e.g. deep inelastic scattering. On the other hand, on the low energy scales of hadron physics the coupling is strong and perturbation theory invalid. The only known nonperturbative and first principles method to compute QCD predictions in this regime is by simulations of lattice gauge theory, whose results are beginning to give a quantitative description of the hadron spectrum.
At a critical temperature T c ∼ 200 MeV, QCD predicts a phase transition to a deconfined plasma of quarks and gluons. At the same temperature the weakly broken chiral symmetry, responsible for the existence of light pions, gets restored. Such temperatures were realized in the early universe, which during its cooling expansion passed through the plasma phase and the quark hadron transition on its way to its present state. We are currently witnessing exciting attempts to recreate this primordial plasma and its subsequent cooling through the phase transition in heavy ion collision experiments at RHIC (BNL) as well as SPS and LHC (CERN). These studies will have a bearing far beyond QCD in the context of early universe physics. Many prominent features of the observable universe, such as the baryon asymmetry or the seeding for structure formation, have been determined primordially in hot plasmas described by non-abelian gauge theories. The QCD plasma is prototypical for those as it is the only one we can hope to produce in laboratory experiments. Moreover, for high densities and low temperatures more exotic non-hadronic phases like a color superconductor have been predicted, and there is a chance for such physics to be realized in the cores of compact stars.
For such applications we need to know how the properties of QCD change under extreme conditions of temperature and density. This entails a determination of the QCD phase diagram and the associated critical properties, a quantitative understanding of the equation of state, the way in which hadron properties get modified as well as the nature of the lightest excitations in non-hadronic phases. Answers to those questions may also help to elucidate mechanisms of confinement and chiral symmetry breaking.
In this review we collect data from lattice simulations of thermal QCD and try to place them into their physics context. While avoiding technical issues as far as possible, our aim is to convey a picture of what we have learned about the QCD transition and plasma, even though in places this picture is still incomplete. When the critical temperature is crossed, light degrees of freedom appear to be released, in accordance with the naive expectation of deconfinement. However, the resulting plasma is not compatible with the perturbative picture of a free parton gas for any temperatures of interest. Rather, the system displays strong residual intercations caused by soft modes in the gauge sector.
Static equilibrium physics is well under control in the continuum limit for the pure gauge theory without quarks, whereas calculations with dynamical quarks are limited to m π > ∼ 300 MeV, and even those are still affected by cut-off effects. Nevertheless, these systematic errors can be reduced step by step in additional calculations to come.
Less advanced and not yet fully quantitative are the calculations in two other areas of interest. Since lattice QCD works in euclidean space-time, it is much harder to obtain a real time description of dynamical processes. This field is in its infancy and we will report on a first step taken in that direction. In the last two years significant progress has been made towards simulating the regime of small baryon densities relevant for heavy ion collisions, which was entirely impossible before then. While this also is a rather young field of research with a long way to go, there now are first results for the (µ, T ) phase diagram and the search for a potential critical endpoint of the deconfinement transition. Before presenting lattice results available today, we give a brief introduction to the theoretical set-up of lattice thermodynamics and a discussion of its current limitations.

QCD Thermodynamics on the Lattice
Detailed introductions to thermal field theory and lattice field theory can be found in [1] and [2], respectively. The QCD grand canonical partition function is given by with the euclidean gauge and fermion actions Its thermodynamic parameters are the temperature T , the volume V and the chemical potential µ for quark number Q (the one for baryon number is µ B = 3µ), while the QCD action depends on the number of quark flavors N f , their masses m f q and the gauge coupling g. In the following we will consider mostly two and three flavors, and always take m u = m d . The case m s = m u,d is then denoted by N f = 3, while N f = 2 + 1 implies m s = m u,d .

The Lattice Formulation
The theory is discretized by introducing a euclidean space-time lattice L 3 ×N t with lattice spacing a, such that volume and temperature are V = (aL) 3 The fermion fields live on the lattice sites x, whereas the gauge fields are represented by link variables U µ (x) ∈ SU(3) connecting the sites. After a suitable discretization of the actions, the Gaussian integral over the quark fields can be performed leading to the lattice partition function where M(µ) = γ µ (D µ + m f − µγ 0 ) is the fermion matrix and S g [U] denotes the gauge action with lattice gauge coupling β = 6/g 2 . For µ = 0, the theory in this form is amenable to a stochastic calculation of expectation values by Monte Carlo methods. Simulations at µ = 0 have only begun to become possible in certain limits, as will be discussed in detail in Section 9. The lattice spacing depends on the gauge coupling, a = a(g), and is set by some physical (zero temperature) quantity like a hadron mass. The physical temperature is then obtained via Eq. (3) as T /m H = 1/(am H N t ). In general, calculations are affected by finite size and cut-off effects, the latter depending on the particular discretization chosen. The task then is to compute observables for various volumes and lattice spacings, and extrapolate to the thermodynamic (V → ∞) and continuum limits (a → 0, N t → ∞), while keeping T and physical parameters fixed.
In practice calculations are severely limited by computational resources, the most expensive part being the evaluation of the fermion determinant. For this reason it has often been set to 1 in the past, corresponding to the quenched approximation which neglects all quark loop contributions. In dynamical simulations at finite temperature, the Wilson and staggered fermion actions are used predominantly, as they are the cheapest to simulate. Another choice with better chiral properties are domain wall fermions, while dynamical overlap fermions so far have not been employed.
Typical lattice sizes of current dynamical simulations are 16 3 ×8 or 32 3 ×4, but for many exploratory studies smaller lattices have to make do. In order to have reasonably small finite size and cut-off effects, the correlation length ξ of a typical hadronic state has to fit comfortably in the box while being much larger than the lattice spacing, Around the deconfinement transition T c ∼ 200 MeV ∼ (1fm) −1 , the temporal size N t = 4 − 8 implies lattice spacings of a ∼ 0.1 − 0.3 fm, allowing for box sizes of aL ∼ 1.5−3 fm. In the chirally broken phase the right inequality then constrains the lightest feasible quark masses, implying that the physical pion cannot be accomodated on such lattices. Typical calculations correspond to m π ∼ 300 MeV. Above T c meson correlation lengths typically scale as ξ ∼ 1/T so that the constraint of Eq. (5) reads Cut-off effects are then controlled by 1/N t independent of temperature. These cut-off effects can be rather severe for certain thermodynamic quantities. The reason is that thermal distribution functions are peaked for modes with momenta p ∼ T of the order of the cut-off, or in lattice units ap ∼ 1/N t . Control over these effects thus warrants the use of improved actions, which are augmented by additional terms irrelevant in the continuum limit and tuned such that, at finite lattice spacing, they subtract cut-off effects to a given order (a n ) in the lattice spacing. In addition, cut-off effects can be smaller on anisotropic lattices, for which a finer spacing in the temporal direction is chosen, a t < a s , corresponding to a larger number of points N t . However, soft gluonic modes scale as ξ ∼ 1/g 2 (T )T , upon which 1/N t → g 2 (T )/N t in Eq. (6), thus requiring large aspect ratios L/N t in both cases, which limits simulations of very large temperatures.

Effective High T Theory: Dimensional Reduction
At large T , when the gauge coupling g(T ) is sufficiently small, a hierarchy between different relevant scales of thermal QCD develops, The lowest non-vanishing bosonic Matsubara mode ∼ 2πT is characteristic for non-interacting particles. The dynamics generates the Debye scale m E ∼ gT , over which color electric fields are screened, and its non-perturbative analogue m M ∼ g 2 T for color magnetic fields [3]. For physics on scales larger than the inverse temperature, |x| ∼ 1/gT ≫ 1/T , this allows an effective theory approach in which the calculations are factorized: integration over the hard modes may be performed perturbatively by expanding in powers of the ratio of scales gT /(2πT ) ∼ g/(2π). This includes all non-zero Matsubara modes, in particular the fermions. It results in a 3d effective theory for modes ∼ gT and softer, and is known as dimensional reduction [4]. Since 4d euclidean time has been integrated over, A 0 now appears as an adjoint scalar, and the effective parameters are functions of the original ones, Discretization and simulation of this reduced problem is easy. Without fermions and one dimension less, much larger volumes and finer lattices can be considered. Moreover, 3d gauge theories are superrenormalizable and the coupling scales linearly with lattice spacing. Hence, very accurate continuum limits can be obtained and systematic errors from simulations are practically eliminated. However, the reduction step entails two approximations: the perturbative computation is limited to a finite order in g and neglects higherdimensional operators, which are suppressed by powers of the scale ratio. The reduction step has been performed up to two-loop order [5] at which parameters have relative accuracy O(g 4 ), while for correlation functions the error is [6] δC/C ∼ O(g 3 ). In the treatment of the electroweak phase transition, this is less than 5% [7], for QCD it depends on the size of the coupling g(T ). As we will discuss in Section 5, it can be non-perturbatively checked that the effective theory accurately describes static correlation lengths at temperatures as low as T > ∼ 2T c , thus allowing for a straightforward treatment of very large temperatures as well as detailed dynamical investigations in the plasma phase.

The (m u,d , m s ) Phase Diagram
On quite general grounds it is expected that at high temperature or baryon density QCD undergoes a transition from hadronic matter to the quarkgluon plasma. However, whether the transition is characterized by a truly singular behavior of the partition function leading to a first or second order phase transition, or whether it merely is a crossover with rapid changes in some observables, crucially depends on the values of the quark masses. At a transition the details of the microscopic interactions are negligible compared to the dominant long range correlations, which are determined by global symmetries. Expectations on the nature of the transition have thus been derived from studies of simpler systems with the same global symmetries as QCD [8], see Figure 1.
The method to locate a transition usually is to look for rapid changes of order parameters like the Polykov loop L (in the pure gauge case) or the chiral condensate (in the chiral limit) and for peaks of their fluctuations, i.e. susceptibilities as e.g. χ L = V (L − L ) 2 . The locations of these changes or peaks define (pseudo-)critical couplings β c , which can be turned into temperatures through the knowledge of a zero temperature quantity like a hadron mass am H (β c ) at the same coupling, T c /m H = 1/(N t am H (β c ). The nature and critical properties of the transition can be obtained from scaling dependencies of various quantities on volume or external fields like the quark mass.
In the SU(3) pure gauge theory, the expected first order transition was  [22] no two-state, Columbia [22] crossover like, JLQCD [20] 1st order like, JLQCD [20] two-state, JLQCD [20] m crit , 3-flavor reweighted [23] m crit , cumulants [19,23]  verified numerically [9] quite a while ago. In the chiral limit of two-flavor QCD the relevant global symmetry group most likely is SU R (2) × SU L (2) ≃ O(4). In analogy to the O(4) σ model, one therefore expects [10] a second order phase transition with O(4) critical properties in the chiral limit. At non-vanishing quark masses the chiral symmetry is broken explicitely and a crossover is expected. All existing lattice calculations [11,12,13,14,15] support a continuous transition. However, O(4) values for the critical exponents could not yet definitively be identified. Alternatively, the anomalous U A (1) could become restored effectively through the disappearance of topologically non-trivial zero modes of the Dirac matrix at or even below the critical temperature rendering the transition first order [10]. This case would be reflected in a degeneracy [16] between the pion and the scalar isovector meson (δ/a 0 ), which does not seem to be realized [17]. For three degenerate massless quarks, again from symmetry reasons, one expects a first order transition. The first order property of the transition should extend from the chiral limit to an end point at some non-zero critical quark mass m c q , where the transition becomes continuous and is expected to be in the Ising universality class [18]. A recent lattice simulation [19] supports this picture even with respect to the universality class. The location of the end point, however, is not a universal property and the current estimates [19,20,21] vary between 190 and 260 MeV for the corresponding pseudoscalar mass, depending on the particular fermion action chosen.
For non-degenerate quarks, i.e. a strange quark heavier than the two light ones, m s = m u,d , the end point becomes an end line separating the first order from the crossover regime, see Figure 1. In the vicinity of the degenerate three-flavor end point the slope of the line can be obtained from a Taylor expansion, m c s = m c + 2(m c − m u,d ). In the currently accessible quark mass range, this relation is consistent with all available lattice data [19,20,22,23], see Figure 1. If one assumes that such a linear extrapolation holds down to the physical light quark mass values, the resulting ratio m c s /m phys u,d ≃ 5 − 10 would put the point of physically realized QCD, m phys s /m phys u,d ≃ 20 (see e.g. [24] for current lattice results), into the crossover regime of the phase diagram.
The location of the phase transition, i.e. its critical temperature, has been known rather precisely for quite a while in the pure SU(3) gauge theory [25,26]. It has since then been confirmed by studies which explored the efficiency of various improved actions [28,29,27], or anisotropic lattices [30,31], in reducing finite lattice spacing effects. The number is most readily given in terms of the dimensionless ratio to the string tension, √ σ ≃ 425MeV, as which was obtained as a weighted average over all available lattice data. For QCD including light dynamical quarks the current results for the critical temperature are summarized in Figure 2. The left panel contains data [32,33,34,35,36] for two flavors, all obtained with improved actions. Discrepancies between Wilson and staggered quarks seen in earlier computations with standard discretizations are greatly reduced. T c is plotted versus the zero-temperature pseudoscalar mass as a measure of the quark mass, both in units of the vector meson mass. The strong decrease of the ratio T c /m V with increasing quark mass is due to the rising vector meson mass which diverges in the infinite quark mass limit, m P S /m V → 1. Close to the chiral limit two conflicting quark mass dependencies should influence the ratio T c /m V : the vector mass should increase linearly with the quark mass, leading to a decrease of the ratio with √ m P S . This is counteracted by the quark mass dependence of T c itself, rising proportional to m 1/βδ q , which amounts to ∼ m 1.1

P S
for O(4) values for the critical exponents. Thus, at small enough quark and pseudoscalar mass, respectively, a linear decrease towards the chiral limit is expected to win. This effect is not yet seen. A scale for T c much less affected by the quark mass is the string tension which is used on the right panel of Figure 2. Indeed, this figure shows that the critical temperature is decreasing with decreasing quark mass. The decrease appears to be compatible with a linear dependence on the pseudoscalar mass as expected for N f = 2. The linear dependence, however, also seems to hold for three degenerate quarks. Note also, that T c is reduced half way down from the quenched value, shown as the band to the very right of the figure, already at a pseudoscalar mass as large as 1500 MeV. It is thus not the chiral dynamics which controls the critical temperature, and one may speculate that a resonance gas picture is more appropriate to describe the thermodynamics close to T c . The overall quark mass dependence for both two and three flavors is rather weak. The difference between N f = 2, 3 is small and amounts to about 20 MeV, independent of the quark mass. Extrapolating the results to the chiral limit, one obtains 2 − flavor QCD : where m ρ has been used to set the scale. Although the agreement between staggered and Wilson quarks is striking, one should keep in mind that the errors are statistical only and do not account for the as yet unknown systematic errors originating from a non-vanishing lattice spacing. One might hope, though, that these are small since improved actions have been used in both simulations.

Equation of State
Energy density ǫ(T ) and pressure p(T ) as a function of temperature are certainly among the most fundamental thermodynamic quantities of the quark gluon plasma governing, e.g., the time evolution of the plasma once being created in heavy ion collisions. At high temperature the relevant partonic degrees of freedom have momenta of order πT which are strongly affected by the UV cut-off introduced through the finite lattice spacing. It is thus important to gain control over the discretization effects and carry out the continuum limit a → 0. Again, this was first accomplished in the pure gluon theory [25], later confirmed in studies utilizing improved actions and/or anisotropic   [32,33,34,35,36]. Right: N f = 2, 3 for improved (p4) staggered quarks, compared to N f = 2 standard (std) simulations [32].
lattices [28,29,31]. As calculations in the ideal gas limit of lattice perturbation theory reveal, the discretization effect is even more pronounced in the presence of dynamical quarks. Since on the other hand the signal for these quantities vanishes proportional to a 4 the use of improved actions, which are designed to reduce the UV cut-off effects, seems mandatory. Moreover, the simulations in the quenched approximation have shown that this improvement helps to extract the continuum limit also in the intermediate T range between 2 and 4 T c .
The results of a computation of the pressure with an improved action in the gauge as well as in the fermion sector [37] are shown in Figure 3. The data have been obtained for N f = 2, 3 with (bare) mass m q /T = 0.4 as well as for N f = 2+1 with a heavier mass m s q /T = 1. For comparison the figure includes the continuum extrapolated quenched result. The figure shows a rapid rise of the pressure in a transition region. The critical temperature as well as the magnitude of p/T 4 reflect the number of degrees of freedom liberated at the transition. In fact, when the pressure is normalized to the appropriate Stefan-Boltzmann limits shown as arrows in the figure, the function p(T /T c )/p SB turns out to be almost flavor independent. It is thus well described by where the ideal gas limit only is modified by an apparently flavor independent function f (T /T c ). At this stage it should be pointed out, though, that the full QCD results have not yet been extrapolated to the continuum limit.  [37]. Right: energy density for N f = 2 improved Wilson quarks on N t = 4 (filled symbols) and N t = 6 (open symbols) [39]. Marks on the right side denote Stefan-Boltzmann limits.
However, the experience gathered perturbatively and in the quenched approximation leads one to expect that the finite a effects will distort these findings by not more than 10%. The most recent results on the energy density are shown in Figure 3 to the right. Note that the dependence on the quark mass appears to be rather weak over a wide range of values. The data has been obtained for N f = 2 improved Wilson quarks [39] on N t = 4, 6 lattices. However, in the particular discretization chosen the improvement does not reduce the infinite temperature finite lattice spacing effects. This is visible in the figure by the big difference in the Stefan-Boltzmann limits computed on the two lattices with different N t values. On the other hand, the cut-off effects are not important close to T c as a comparison between N t = 4 and 6 results reveals. This is expected since in this temperature regime the correlation length is large and infrared modes dominate. Thus, a comparison of the different fermion discretizations is meaningful. In fact, improved staggered data are available [40] and lead to a consistent estimate for the critical energy density of ǫ(T c ) ≃ (6±2)T 4 c , which also agrees with the value originating from an earlier simulation using the standard staggered action [41].
In all cases, up to the highest temperature investigated, pressure and energy deviate substantially from the ideal gas behavior. The deviation is too big to be reproduced in ordinary high temperature perturbation theory which converges badly in this temperature range [42]. A more refined tech- nique as hard thermal loop resummation [38] gets closer to the lattice results. The question of the approach to the ideal gas limit at higher temperatures can also be addressed by dimensional reduction, which splits the thermodynamic functions into perturbatively calculable parts and a remainder to be simulated in the effective theory Eq. (8) [43]. Unfortunately, this requires determination of an integration constant e 0 by a four-loop calculation, which has not been completed yet. However, results [44] treating e 0 like a free parameter and matching it to the 4d simulations are shown in Figure 4. If the eventual value for e 0 lies in the right range, the whole picture makes sense and explains how ideal gas values are only attained at asymptotic temperatures.

Screening Masses
Essentially all static equilibrium properties of a thermal quantum field theory are encoded in its equal time correlation functions. Even though these quantities may not be directly accessible in heavy ion collision experiments, their theoretical knowledge provides us with the relevant dynamical length scales in the plasma, from which conclusions about the acting degrees of freedom and their physical effects may be drawn. The connected spatial correlation functions of gauge-invariant, local operators A(x), fall off exponentially with distance. The "screening masses" M are the eigenvalues of the spacewise transfer matrix of the corresponding lattice field theory, and classified by its symmetries. Because of the shortening of the euclidean time direction at T > 0, the continuum rotation symmetry of the hypertorus orthogonal to the correlation direction is broken down from O (3) to O(2) × Z(2), and its appropriate subgroup for the lattice theory is D 4 h . The irreducible representations and the classification of operators have been worked out for pure gauge theory [45,46] as well as for staggered quarks [47]. Physically, the screening masses correspond to the inverse length scale over which the equilibrated medium is sensitive to the insertion of a static source carrying the quantum numbers of A. Beyond 1/M, the source is screened and the plasma appears undisturbed. Technically, the computation of these quantities is equivalent to spectrum calculations at zero temperature.  Figure 6) shows that the pure gauge results are close to continuum physics. In the range 0.8T c < T < T c , the lowest scalar screening mass is observed to be roughly 20% lower than the lightest scalar glueball at zero temperature, M(T )/M G (T = 0) ∼ 0.8. At T c a sharp dip is observed, after which the screening masses appear to be proportional to T . Screening states with larger masses show the same qualitative behavior above T c .

Hadronic Screening Masses
Apart from operators A made from pure glue in Eq. (11), also spatial correlations of meson operators have been investigated, both in the quenched approximation as with various numbers of dynamical fermions [49,50,51,52,53,54,55,56,57]. The picture which has emerged so far is illustrated by some of the available data shown in Figure 5 (right). Below T c , the screening masses have shown neither a marked temperature dependence nor, correspondingly, a drastic difference to zero temperature masses. At temperatures above T c spatial (as well as temporal) correlation functions reflect the restoration of the chiral SU L (N f ) × SU R (N f ) symmetry. In particular, the vector and axial vector channel become degenerate independent  [48]. Right: Mesonic screening masses in the quenched chiral limit. Below T c , M ρ /T c [51,52] is plotted, the line denotes the T = 0 value. Open Wilson symbols denote anisotropic lattices [54]. Staggered pion data are extrapolated to a = 0 [56], Wilson and staggered rho are from N t = 8 [51] and 16 [52] lattices. The free quark limit has not been corrected for finite volume effects.
of the discretization and of the number of dynamical flavors being simulated. Moreover, the pion ceases to be a Goldstone boson and acquires a screening mass. A degeneracy is further observed within errors in the pseudoscalar and isoscalar scalar (σ/f 0 ) channel although the latter, for technical reasons, is difficult to access on the lattice. Nevertheless, screening masses obtained from fits to correlation functions [57] and susceptibilities (at finite lattice spacing) [58] lead to a consistent picture. At T c the anomalous U A (1) most likely is not yet restored effectively [17], see also [53,57,58,59], as has been mentioned already in Section 3.
At high temperature, the screening masses are expected to approach the free quark propagation limit, M → 2πT . In fact, already at temperatures of about 1.5 T c , the results for the vector channel are not far from this value. If the finite volume effects for free quark propagation are taken into account, the deviations amount to about 15 % and decrease only slowly with temperature. Moreover, at temperatures > ∼ 1.5T c , in the Wilson discretization pion and rho become nearly degenerate already at finite lattice spacing. Earlier staggered simulations had found discrepancies here. However, a recent paper [56] reports that also in the staggered discretization a degeneracy of pseudo- scalar and vector is reached in the continuum limit. Quenching effects are found to be below 5% for T > T c [60].
In the plasma phase, screening masses can also be investigated within the dimensionally reduced theory, cf. Section 2.2. In this framework they correspond to the spectrum of the transfer matrix for the 3d theory. The associated Hamiltonian respects SO(2) planar rotations, two-dimensional parity P , A 0 -reflections R, and the symmetry is again SO(2) × Z(2) × Z(2) = O(2) × Z(2). Remember however that in this setup one is interested in soft modes, while the Matsubara frequency ∼ 2πT represents the UV cut-off for the effective theory. Figure 6 displays numerical results for the lowest screening masses for the SU(2) and SU(3) pure gauge theories, obtained through simulations of both the 3d effective theories and the full 4d theories, in various quantum number channels. Even at a temperature as low as T ∼ 2T c , good agreement is observed except in one case. Note that the error bars are only statistical, while both calculations have systematic errors: the effective theory has truncation errors, while the 4d simulations are not infinite volume and continuum extrapolated. Estimating such effects at 10-15%, we conclude that within current accuracy dimensional reduction gives a quantitative description of screening lengths for temperatures T > ∼ 2T c . The same conclusion is reached with other observables like the spatial string tension [64], static potentials [65] or gauge fixed propagators (cf. Section 5.3). For larger temperatures the M/T get lower and the spectrum denser, as shown in Figure 6 (right) for a few states. A large number of such screening masses in many quantum number channels for various numbers of flavors has been computed in [62].

How Perturbative is the Plasma?
The approach of dimensional reduction is particularly valuable in disentangling contributions from different degrees of freedom, thanks to accurate mixing analyses, as well as for treating larger temperatures T ≫ T c which cannot be reached by 4d lattices. This can be used to inspect to what extent the plasma behaves perturbatively. Figure 6 then teaches us that for any reasonable temperature the largest correlation length of gauge-invariant operators belongs to the A 0 degrees of freedom and not to the A i , in contrast to the naive parametric ordering, Eq. (7). Only asymptotically is the perturbative ordering restored.
Of particular phenomenological relevance for the QCD plasma is the Debye mass, whose inverse gives the length scale over which colour-electric flux is screened. Its expansion in powers of coupling constants reads with the leading order perturbative result m LO D = (N/3 + N f /6) 1/2 gT . At next-to-leading order ∼ g 2 T , only a logarithm can be extracted [66], whereas the coefficient c N is entirely non-perturbative.
A gauge invariant definition of the Debye mass has been suggested in [67], according to which it corresponds to the mass of the lightest gauge invariant screening state odd under euclidean time reflections. In the 3d effective theory the latter is replaced by the scalar reflection R, and the Debye mass thus corresponds to the 0 −+ − ground state ∼ TrA 0 F 12 . The coefficient c N can be measured separately from the exponential decay of a Wilson line in a 3d pure gauge theory [67], which has been performed in [68]. This allows to disentangle the contributions from different scales, as shown in Table 1. The O(g 3 T ) corrections are less than 30% even at temperatures as low as T = 2T c , and they disappear entirely for asymptotically large temperatures. On the other hand, the dominant scale is again ∼ g 2 T for all Table 1: Contributions of the first two, the third and higher order terms in Eq. (12) to the total m D [62].  (12) temperatures of interest, in contrast to the naive expectation ∼ gT . Only at asymptotically high temperatures is the perturbative picture restored. These findings are compatible with and explain the behavior of the pressure density, which attains its ideal gas limit only at asymptotic temperatures. There the effect is less pronounced, since the pressure is dominated by hard modes and soft modes provide corrections, whereas the screening masses discussed here couple exclusively to soft modes.

Gluon Screening Masses
With energy and pressure densities as well as the quark number susceptibilities rising rapidly across the transition, it seems clear that some light constituents are liberated in the deconfined phase. A longstanding question concerns the nature of these degrees of freedom and their properties. Moreover, spatial correlators of singlet operators as discussed in the previous subsections, and in particular a Debye mass defined by one of them, are not directly related to phenomena like J/ψ-suppression. These are caused by color charged intermediate states, i.e. some kind of constituent. In perturbation theory one therefore defines the Debye mass in analogy to QED as the pole in the A 0 -propagator [66], with a perturbative expansion as in Eq. (12), which can be proved to be gauge invariant order by order [69]. On the other hand, at moderate temperatures T > T c interactions are still strong, and a purely perturbative parton picture is surely not appropriate. Similarly, the propagator of magnetic fields A i is expected to develop a magnetic mass on the scale ∼ g 2 T , which however even at asymptotic temperatures is entirely non-perturbative [3]. Unfortunately, lattice studies of gluon propagators [70] are hampered by several problems. It is difficult to fix a gauge uniquely and avoid the problem of Gribov copies [71]. Moreover, most complete gauge fixings, like Lan-  (2) gauge theory. The line represents a fit to 4d simulations [73], the symbols show simulations of the dim. red. theory with different matching procedures. Triangles correspond to two-loop perturbative matching. From [74].
dau gauge, violate the positivity of the transfer matrix, thus obstructing a quantum mechanical interpretation of the results. However, recently it was shown analytically [72] that gluon correlators can be made gauge invariant by dressing them with appropriate (gauge fixing) functionals. These correlators fall off exponentially with eigenvalues of the Kogut-Susskind Hamiltonian in the presence of sources [72], thus probing the energy of a gluon coupled to sources. The spectrum is gauge invariant, while matrix elements depend on the particular functionals chosen. This is true for all gauge fixings which are local in time (to preserve the transfer matrix), such as the Coulomb gauge. Numerically, in detail only SU(2) pure gauge theory has been investigated so far. Figure 7 shows results obtained for the electric screening mass obtained both in 4d simulations as well as in simulations of the dimensionally reduced theory, where a variety of gauges was found to agree with the Coulomb gauge results. Like for the hadronic screening masses, with two-loop matching good agreement between full and effective theory for temperatures down to T ≈ 2T c is observed. Equally consistent is the finding that the leading ∼ gT contribution to the electric mass is subdominant for temperatures up to O(10 4 T c ) [73].
Extracting the magnetic mass appears to be more difficult technically. In Landau type gauges, where positivity is violated, severe finite volume effects are observed [73,74]. In Coulomb type gauges the spatial non-locality of the fixing functionals leads to an exclusive projection on torelonic states [72]. However, by introducing explicit local fields to mimic static sources and extrapolating to their infinite mass limit, the magnetic mass has been calculated in the high temperature limit, where the dimensionally reduced model is just the 3d SU(2) pure gauge theory, and the result is [75] m M = 0.36(2)g 2 T.
It is interesting to note that this compares well with various resummation methods leading to self-consistent gap equations, which at one loop yield pole masses in the range of m M ∼ 0.25−0.38g 2 T [76], while a two-loop calculation gives m M = 0.34g 2 T [77]. Note that these partonic screening masses are much lower than those of the singlet operators, and hence correspond to larger correlation lengths. This should play an important role for fragmentation phenomena in the plasma. For a constituent picture of the plasma to be appropriate, it should be possible to interpret the singlet operators as multigluon states [45]. For high T within the 3d theory, this means [78] that the masses of the lowest dimensional operators should be approximated by appropriate multiples of m E , m M . This appears to work up to about 10% for temperatures T > ∼ 100T c [78,79], but accuracy decreases with T . E.g. the two masses shown in Figure 5 (left) seem to slowly approach their multigluon ratio 3/2 from above [48].

Free Energy of Static Quarks
A different source of information on gluonic excitations is provided by heavy quark free energies, in particular of a quark-antiquark pair. Such systems are also of interest for the physics of heavy quarkonia in the medium, in particular for the question of J/ψ-suppression [80]. The QQ free energy is defined [81] by the partition function of the thermal heat bath containing a static quark and antiquark source at separation R, where trL( At zero temperature the free energies reduce to the heavy quark potentials, and at non-zero temperature they also exhibit modifications of the potentials by Boltzmann weighted thermal excitations. In the pure gluon theory, at temperatures below T c and large distances R, the free energy rises linearly. The coefficient of the linear term, the string tension, decreases with increasing temperature and vanishes above T c [84]. In the presence of dynamical quarks the color charges of the heavy quarks are screened also below T c and one observes [85,32,86] the expected string breaking, Figure 8. The distances where the free energies become flat in R range from 1.5 to 1 fm, decreasing with temperature, even at (bare) quark masses as large as m q /T = 0.4. Note that the deviations from the zero temperature quenched potential set in already at distances of O(0.5fm) for temperatures > ∼ 0.75T c . Normalizing the free energy to the short distance zero temperature potential, its large R asymptotic value rapidly decreases with T .
In the deconfined phase the free energy shows the behavior of an exponentially screened potential, and is well described by where c(T ), d, µ(T ) are fit parameters. In perturbation theory, the leading term originates from two-gluon exchange and predicts d = 2 and exponential decay with twice the Debye mass [82]. Lattice investigations [84,73] have shown that this simple behavior does not apply in the temperature range explored. Rather, fits [84] to Eq. (16) favored d ≃ 3/2 and found screening masses µ(T ) to be compatible with the lowest color singlet 0 ++ + screening mass shown in Figures 5, 6. This is not surprising: the Polyakov loop is a gauge invariant operator, and since it is an exponential of gauge fields it couples to all J P C sectors. Consequently its correlator decays with the lightest screening mass of the spectrum.
In order to address Debye screening of a static quark-antiquark system, one needs to access the singlet free energy channel, which in the deconfined phase is expected to decay exponentially with the electric mass m E [83]. Constructing thermal correlators for the different channels separately in general requires gauge fixing. Recently it has been shown that this can be achieved in a manner respecting the transfer matrix, and such correlators decay with a gauge invariant spectrum [87]. The octet channel was found to be repulsive at short distances, as predicted by perturbation theory [88,82], while at large distances it approaches the confining singlet potential from above. First studies of the thermal singlet channel have been reported in [89].

Susceptibilities
Susceptibilities of quark number and other quantities offer an elegant and simple way to study effects of small finite baryon densities on the equation of state within simulations performed at zero chemical potential [90]. Moreover, some of them are directly related to event-to-event fluctuations of particle production in heavy ion collisions, and help to providing signals to differentiate between the phases in order to detect the quark gluon plasma [91].
susceptibilities in the two flavor case, χ sing = 2χ uu +2χ ud , χ trip = 2χ uu −2χ ud . Numerical studies with dynamical staggered quarks show a rapid jump at the critical temperature [92,53]. A more recent study [93] for 2+1 flavors is shown in Figure 9, where the strange quark mass was fixed by tuning m ηss /m φ to its physical value, i.e. m M S s (2GeV) ≈ 95 MeV [24]. In addition it contains the diagonal strange quark number susceptibility, which at large temperatures should become equal to half the triplet susceptibility. The use of improved actions (one loop Symanzik for gauge fields and Asqtad for fermions) leads to good apparent scaling behavior, suggesting deviations from the free quark gas value at 2T c of not more than 10 % in the continuum limit.
Other studies [60,56] considered mostly the deconfined phase, T ≥ 1.5T c . They find that in this regime quenching only affects susceptibilities by about 5%, while cut-off effects for unimproved actions play an important role. Based on these observations, quenched calculations using two different fermion discretizations were used to obtain continuum extrapolated results [56], as shown in Table 2. These authors also study susceptibilites for baryon number, χ B = (4χ 3 + χ s + 4χ ud + 4χ us ), and electric charge, χ Q = (10χ 3 + χ s + χ ud − 2χ us ), where χ 3 is χ trip /2. Flavor off-diagonal susceptibilities are found to be consistent with zero, which is in contrast with Table 2: Results for the continuum limit of quark number susceptibilities in quenched QCD, extrapolated from staggered quark simulations with valence mass m s /T c = 1. From [56]. (1) 1.09 (3) predictions from a three-loop perturbative calculation [94]. It is also inconsistent with the results of [93], which reports statistically significant off-diagonal elements of the order of 10 −2 at T ∼ 1.5T c . While this issue needs clarification, off-diagonal values in any case seem to be small compared to the diagonal ones except very close to T c [93]. For the isovector susceptibility χ 3 numerical results are compatible. At T = 3T c , χ 3 from Table 2 falls short of its free fermion gas value χ 3 F F = T 2 by about 10%. For such temperatures and larger, it is consistent with HTLresummations, which predict χ 3 /χ 3 F F ≈ 0.90 − 0.94. This observation fits to the behavior of the pressure and the discussion of screening masses at comparable temperatures, with soft gluon modes accounting for remnant interactions. It is interesting to note that for the isovector chemical potential µ 3 = µ u − µ d the fermion determinant is positive, thus it can be simulated (cf. Section 9) and numerical results can be compared to analytic predictions [95].
In a similar vein, responses of the lowest pseudoscalar screening mass to scalar and vector chemical potentials have been investigated in a dynamical simulation with two light staggered flavors [96]. Because of the symmetry of the grand canonical partition function under combined euclidean time and µ-reflections, terms odd in (µ/T ) vanish for screening masses, which is also observed numerically. The result for the leading quadratic coefficient is shown in Figure 9, where µ S = µ u = µ d , and µ V = µ u = −µ d . Again one observes a pronounced jump across the deconfinement transition. Below T c , the response to a small scalar chemical potential is essentially zero, consistent with the fact that pions are still Goldstone bosons associated with chiral symmetry breaking in this regime. This ceases to be the case above T c , where a large response is observed. The behavior with isovector chemical potential is quite different, and consistent with expectations based on the phase structure conjectured in [95].

Dynamic Properties
Lattice simulations inevitably are carried out in euclidean space-time. In order to make predictions for real-time processes, in general an analytic continuation to Minkowski space has to be performed. At zero temperature, this is trivial for certain quantities such as masses of stable states or certain matrix elements. At finite temperature the situation is more complicated [1], since e.g. the momentum space propagator is defined only at discrete Matsubara frequencies iω n . Furthermore, because Lorentz invariance is lost due to the presence of the heat bath, spatial and temporal correlation functions are different in general, and the static results as discussed in Section 5 can not immediately be used in the interpretation of experimental results for dynamical processes.
The full information about plasma excitations, the presence of genuine particle poles, resonances, their location and widths etc., is contained in the spectral density σ H (p 0 , p) in a given channel with quantum numbers H. It is related to temporal correlation functions by Temporal correlations have been used in order to gain information about the spectral density. However, since the extension of the system in the temporal direction is limited by the inverse temperature, the isolation of a ground state dominating the correlation function at large distances is difficult. To improve the reliability of fits, extended (smeared) operators of various kinds have been used, sometimes combined with anisotropic lattices [50,54]. With smearing one hopes to increase the projection onto the ground state, whose contribution would then dominate already at small temporal distances. Anisotropic lattices increase the number of data points in the temporal direction, to stabilize more sophisticated fit ansätze including ground and first excited states. However, one still has to rely on fit ansätze and, if possible, use the quality of the fit to distinguish between various models. Moreover, modeling the operators introduces bias, and the results have thus turned out to depend on the method. However, at least qualitatively the results are in agreement with a significant two-quark cut contribution, while in addition the behavior of wave functions was interpreted to indicate meta-stable bound states in the plasma phase [54].
Clearly, it would be much more preferable to obtain the spectral density σ H directly. However, since the lattice provides only a limited set of noisy data points at discrete values of τ , its numerical extraction from correlation functions is an ill-posed problem. Here, progress has been made recently [97] by applying the maximum entropy method (MEM) [98], which attempts to construct the most probable spectral density given the data and taking into account prior knowledge, such as the positivity of σ H and its perturbative behavior at large energies. The method has been successfully tested and applied at zero temperature [97,99], where the contributions of several excited states to the spectral density could be established. At finite temperature one still is hampered by the limited temporal extent of the lattice. Moreover, it is important to check whether one is able to reproduce the two-quark cut in the free case [100], because this is expected to contribute significantly also in the interacting case at high temperature. First results for the vector channel at small quark masses [101] are shown in Figure 10. The main feature of the data is that the rho peaks at low temperature disappear and develop into broad "resonances" above T c whose locations move proportional to the temperature. The vector spectral density σ V (p 0 , 0) is immediately related to the thermal cross section for the production of dilepton pairs at vanishing momentum [102], This is shown in Figure 10 (right). The "resonance" like enhancement of σ V results in the enhancement of the dilepton rate over the perturbative tree level (Born) rate [103] for energies p 0 /T ∈ [4,8]. In the low energy range the spectral density decreases rapidly, in contrast to perturbative calculations which have σ V diverging in this limit [102]. Fortunately, from inspecting the correlation function directly [101] it is already clear that the spectral density vanishes ∼ p α 0 with some power α in the p 0 → 0 limit. However, gaining detailed control over the low energy behavior of the spectral density is a challenge at this point and will become even more difficult at higher temperatures. On the other hand, the zero energy limit of spectral functions From [101].
is related to transport coefficients [104], and it would be extremely interesting to obtain results for these quantities [105] as a step into first principles investigations of non-equilibrium properties of QCD.
9 Lattice QCD at Finite Density QCD at finite baryon density plays a role in two rather different regimes.
In heavy ion collisions the initial state has a small non-zero baryon number, corresponding to baryon chemical potentials of order ∼ O(50 MeV) [106], and any subsequent plasma state is a state of high temperature and low density. On other hand, the core of neutron stars is composed of cold and very dense nuclear matter. These two situations correspond to the regions close to the axes in the tentative QCD phase diagram Figure 11. While in the latter case a rich phase structure including a color superconducting state has been conjectured [108], understanding the former regime is particularly pressing in view of current heavy ion collison experiments. Of special importance is a prediction of the location of the tricritical point, in which the first order deconfinement transition line is expected to terminate. Simulations of lattice QCD have so far failed to be a viable tool for finite density because of the so-called "sign problem". For µ = 0 the relation γ 5 Mγ 5 = M † guarantees positivity of the fermion determinant, det M(0) ≥ 0. However, for SU(3) and µ = 0 this relation does not hold and det M is Figure 11: Conjectured phase diagram for QCD with two massless and one heavy flavour. From [107].
complex. This prohibits Monte Carlo importance sampling, which interprets the measure as a probability factor and thus requires it to be positive. Note that all physical quantities like Z or observables O are real positive, i.e. all imaginary parts cancel out of the path integral and the problem lies merely in evaluating it.
This problem is generic for many fermionic models, and for a certain class of them it can actually be solved by the use of cluster algorithms. These subdivide the ensemble into parts within which all contributions add up positively. The model closest to QCD for which this succeeded is the three state Potts model [109]. Unfortunately, these algorithms are not directly applicable to QCD, and we refer to [110] for an overview.
We will likewise not discuss the considerable amount of work done on certain modifications of QCD for which the sign problem is manageable or altogether absent, such as QCD in the static quark limit, two-color QCD or QCD at finite isospin. These topics are covered in recent reviews [111]. It has also been argued that for very high densities an effective theory free of the sign problem can be derived [112], but no numerical results exist as of yet.
Here we wish to focus on some recent approaches which made significant progress towards simulating QCD with realistic parameter values at small baryon densities. It is important to note that all of them circumvent the sign problem rather than solving it. The partition function in its first form in Eq. (2) suggests that for smallμ = µ/T the problem is close to the trivial case, and in particular lends itself to a Taylor expansion. The following approaches in one way or another all require small enoughμ to succeed.
As we shall see, however, this limitation still allows to obtain results up to µ B ∼ 500 MeV, thus fully including the parameter region relevant for heavy ion collisions. In the following we limit our discussion to mapping out the deconfinement transition line.

Reweighting: The Glasgow Method
A mathematical identity can be used to reweight from configurations generated at zero density, where there is no sign problem, to a corresponding one at finite density. The method is based on rewriting the partition function as where the determinant ratio is now treated like an observable while the integration measure is defined at µ = 0, and hence positive. This expression is exact. However, in a Monte Carlo integration with finite statistics two problems arise. Firstly, the complexity of the numerator leads to oscillations and large cancellations. The reweighting factor corresponds to a ratio of two partition functions with different actions, and thus decays exponentially with the difference between their free energies, which is an extensive quantity, Z = exp −∆F/T ∼ exp (−const.V ). The statistics required for a given accuracy thus grows exponentially with volume, rendering extrapolations to the thermodynamic limit extremely difficult. Secondly, Monte Carlo approximates the partition function by the contributions of the most likely configurations obtained by probability sampling. Since the extrema of the action with µ = 0 are shifted from those of the entire integrand including the determinant ratio, the probability distribution of the generated ensemble is changed. This should be of little consequence for small µ and large statistics, but with growing µ the overlap between reweighted and full ensemble deteriorates. The effect does not show up in the statistical error, and the difficutly is in knowing when this becomes problematic leading to incorrect results. For example, at T = 0 the transition to nuclear matter is found at an unphysical value for µ c . A detailed discussion with early results and references can be found in [113].

Multiparameter Reweighting
Only recently it was realized that the ensemble overlap can be significantly improved by a multidimensional generalization [114] of the Glasgow method.
In addition to µ one may also reweight in the lattice gauge coupling β, so that the partition function is now written as This is crucial for finding critical behavior. While reweighting only in µ would inevitably generate an ensemble away from criticality, the second reweighting parameter can be used to keep it fluctuating between the phases, as the unshifted ensemble certainly would. This method was compared with the Glasgow method in four flavor QCD at imaginary chemical potential [114], for which the determinant is positive and a full ensemble is also available (cf. Section 9.4). Multiparameter reweighting is found to agree significantly better with the full ensemble up to some critical imaginary chemical potential, beyond which it breaks down [115].
The method was then used in [116] to map out the deconfinement line, using Lee Yang zeros [117] as observables. These are points of vanishing partition function, marking the singularities in the free energy. Singular phase transitions only appear in the thermodynamic limit, while on finite volumes the free energy is regular. In the latter case the Lee Yang zeroes are displaced to complex parameter values. As the volume is increased, they approach their thermodynamic limit values on the real axis if there is a genuine phase transition, whereas for a smooth crossover the zeroes settle somewhere in the complex plane. In the case of a first order transition the large volume behavior is consistent with a β c (V ) = β c (∞) + ζ/V behavior. The simulations reported in [116] were carried out for 2+1 flavors of staggered fermions, the result is shown in Figure 12. It constitutes the first numerical prediction of the (µ, T ) phase diagram including a critical point. However, larger volumes are required to check for thermodynamical behavior and confirm the results. Recently attempts have been made to further improve the overlap by splitting the reweighting factor in a product of many factors, and preliminary results on small lattices have been reported in [118].

Taylor Expanded Reweighting
Since for heavy ion collisions one is interested in rather small chemical potentials of a few ten MeV, one may use the smallness ofμ in an approximation of multiparameter reweighting in order to make simulations on larger volumes feasible [119]. Consider a Taylor expansion of the reweighting factor R as a power series inμ = µ/T , and similarly for any operator O. Expectation values are then given by For sufficiently small chemical potentials, this series should converge quickly and calculation of the leading coefficients should give a good approximation to the full answer. The benefit is that the derivatives of the determinant can be expressed through local operators in the fermion fields. These are much cheaper to compute than the whole determinant, which requires a non-local operation. In particular, keeping only the leading term, the cost is equivalent to that of computing susceptibilities (cf. Section 7). Even though the expanded phase of the determinant still has oscillations growing with volume andμ, the saving in computer time can be used to gather higher statistics for larger volumes.
This strategy has been applied in [119] to compute the location of the deconfinement line for a two-flavor theory with p4-improved staggered fermions and a pion mass of about m π ≈ 600 MeV. The method to locate a phase transition is to look for the peak χ max = χ(µ c , β c ) of the susceptibilities where the chiral condensate and the Polyakov loop were used for O in practice. This peak implicitly defines a critical coupling β c (μ), which must be even in µ because of the symmetry of the partition function [115]. This can be employed to calculate the curvature of the critical temperature at µ = 0, Indeed a lattice with twice the spatial size than that in [116] was employed without major problems. Of course, in order to judge the convergence properties of the series, several consecutive coefficients are needed.

Imaginary Chemical Potential
For imaginary chemical potential, the fermion determinant is positive and simulations are as straightforward as for µ = 0. It is therefore natural to ask whether such simulations can be exploited to learn something about real µ. Let µ R , µ I ∈ IR denote the real and imaginary parts of a complex µ = µ R +iµ I . In the presence of a complex chemical potential, a Z(3) transformation of the fermion fields is equivalent to a shift in µ I , causing the partition function to be periodic in µ I with period 2πT /3 [120], Hence, onceμ I exceeds some critical valueμ c I , a phase transition to a nontrivial Z(3) sector occurs, and because of the symmetry Eq. (24) this transition is periodically repeated at the exact critical valuesμ c I = 2π(k + 1/2)/3. It was conjectured [120] and verified numerically [115,121] that these transitions are of first order in the deconfined phase and continuous in the confined phase. On the other hand, for purely real µ these transitions are absent.
Through its periodicity, the grand canonical partition function at imaginary chemical potential is the Fourier transform of the canonical partition function at fixed quark number Q [120,122], It was thus suggested to simulate Z at imaginary chemical potential and perform the Fourier transformation numerically [123]. This has been carried out in the 2d Hubbard model at large T and small Q, but not yet in QCD. Unfortunately, this approach only postpones the sign problem from the Monte Carlo integration to the Fourier transform. As Q is getting large approaching the thermodynamic limit at fixed baryon density, it will again cause rapid oscillations and the associated cancellations.

Analytic Continuation
Ifμ is small another line of attack is possible [124]. Since there are no massless modes in the theroy, all observables are analytic functions ofμ everywhere but on the critical lines of phase transitions. After computing expectation values with Z(iµ I , T, V ), one may therefore fit them by truncated Taylor In cases where this is possible to satisfactory accuracy, analytic continuation to real µ is trivial. There are two advantages to such a procedure. First, there is no sign problem and hence no restriction on the volumes that can be simulated. Second, the simulation gives the value of the whole expectation value rather than just a single Taylor coefficient. While the expansion is necessary in order to analytically continue, fitting to various orders allows control over the convergence properties. However, a check of the convergence and hence controlled continuation are limited toμ I <μ c I = π/3, which marks the first Z(3) transition discussed above.
Non-perturbative evidence for the viability of this approach has been given for screening masses in the deconfined phase [125], which can be simulated sucessfully for imaginary and real µ in the framework of dimensionally reduced QCD [62].
It was shown in [115] that the same approach can be extended to study the critical line of the deconfinement transition itself. Again the transition was defined by the peak χ max = χ(µ c , β c ) of plaquette, chiral condensate or Polyakov loop susceptibilities. On finite volumes these are analytic functions over the whole parameter space of the theory, non-analyticities associated with phase transitions develop only in the thermodynamic limit. The implicit function theorem then implies that also β c (μ) is analytic for allμ, which therefore has a Taylor expansion As the thermodynamic limit is approached, β c (μ) approaches the unique critical line of phase transitions, while remaining pseudo-critical in a crossover regime. In principle the nature of the line can also be determined from its finite volume scaling.
In [115] a two staggered flavor QCD was simulated, and a striking finding is that for allμ I <μ c I = π/3 the critical coupling is indeed well described by a fit quadratic inμ I , with negligible effect of terms O(μ 4 I ). The same conclusion is reported from an imaginary µ simulation of a four flavor theory in [121]. In principle, of course, there is nothing that guarantees such good behavior, which would need to be checked for every change of parameters. It is therefore amusing to note that both in the 3d Gross-Neveu model as well as in some random matrix model with symmetries similar to QCD, the leading quadratic term very well reproduces the exact soultion for the critical line over wide ranges of parameter space [121].

The (µ, T ) Phase Diagram
In Figure 12 we assemble the various phase diagrams published so far in one plot. Some caution shoud be used when interpreting it, as the curves correspond to different parameter values and volumes. The only thing equal in all four simulations is the lattice spacing which at a ≈ 0.3 fm is large enough to cause cut-off effects. Furthermore, the error bars shown on the data points should not be directly compared with the smaller error bands corresponding to fit results.  Table 3.
Nevertheless, the curves are only shown as far as the validity of the quadratic approximation has been checked explicitly in the imaginary µ approach. Within this range it was found in [119] that doubling the quark mass did not affect the line within the current errors, so that comparing different quark mass values makes sense. Moreover, the four flavor calculation at imaginary µ in [121] is entirely consistent with the (unexpanded) reweighting calculation performed at exactly the same parameter values [114]. Table 3: Simulation parameters for the results shown in Figure 12.

Method
N f m q largest lattice reweighting [116] 2+1 am = 0.025, m s = 8m u , (m π ≈ 300MeV ) 8 3 × 4 rew. +Taylor [119] 2 am = 0.1, (m π ≈ 600) MeV 16 3 × 4 imag. µ [115] 2 am = 0.025, (m π ≈ 300) MeV Finally, within the Taylor expanded reweighting it was checked that simulating real and imaginary chemical potential indeed produces the same leading coefficient [126]. These cross-checks are extremely important, since all methods introduce some systematic error growing with µ/T . However, it seems safe to conclude that they all give compatible and reliable results for the pseudo-critical temperature at least up to µ B ≈ 500 MeV. As expected on physical grounds, T c (µ) decreases faster when more light quarks are present. An important physics conclusion then is the flatness of the critical line, which implies that the critical temperature relevant for RHIC is essentially the same as at zero density. It also implies that T c (µ B → 1GeV) has to drop rather dramatically if our qualitative picture of the phase diagram is at all correct. First investigations of the equation of state at finite µ have found negligible effects at the few percent level for RHIC densities [127,119].

The Critical Point, Quark Mass Dependence
The only simulation that has so far produced a critical point is the one employing reweighting on lattices of rather limited size [116]. On the other hand, experience at zero density indicates that rather large volumes are required in order to reach the asymptotic finite volume scaling required for an unambiguous determination of the order of a phase transition [19]. Such systematic effects may well be larger than the statistical error box indicated in the figure, and it is hence important to have another independent calculation. Moreover, while the location of the (pseudo-)critical line appears to be only very weakly quark mass dependent, the location of the critical point strongly depends on them: for three chiral quarks, the first order line joins the temperature axis, while for the 2+1 flavor case depicted in the figure it is quite far out.
The task of locating the critical point as a function of parameters is equivalent to mapping out the critical surface separating first order transitions from crossover in Figure 13, which is simply the extension of the phase diagram Figure 1 to finite densities. Probably the most economical way of doing this is to find the projections of this surface onto the (T, µ) and various (m q , µ)-planes. Just like the former, the (m q , µ)-plane then shows a critical line separating first order from crossover. For this line, m c q (µ), the same analyticity and symmetry considerations apply as discussed for T c (µ), i.e. it has a Taylor series even in µ. Within the Taylor expanded reweighting approach, a first result for the leading coefficient in the case of three degenerate flavors of p4-improved fermions (m π ≈ 200MeV) has been reported in [23] as T ∂m c q /∂(µ 2 ) = 0.21 (6). Moreover, as argued in [128], when the quark mass is varied, the 3 flavor critical endpoint will trace out a line T c (µ c , m q ), which has an analytic continuation to imaginary µ. This in turn allows to determine it from imaginary chemical potential simulations.

Summary
In this article we have attempted to summarize the current understanding of QCD at non-vanishing temperature as it arises from numerical studies on the lattice. While many of the results in the pure gauge theory are available in the continuum limit, simulations with dynamical fermions still suffer from systematic errors. These are mainly due to finite lattice spacing as well as quark masses which don't meet their physical values yet. Nevertheless, for critical temperatures and the equation of state quantitative results could be given, which are expected to be correct at the 10 -15 % level. Apart from the basic thermodynamic quantities, existing results provide us with a detailed picture of how static quantities change through the deconfinement transition up to a few T c . Combinations of perturbative calculations and numerical simulations have produced insight into the regime of very high temperatures as well as the dynamics and mixing of modes. Altogether this led to a quantitative understanding of the relevant static length scales in the plasma, as well as tests of the applicability and breakdown of thermal perturbation theory.
The naive picture of the deconfined phase as a weakly interacting parton gas is not supported. For temperatures relevant to heavy ion collisions, the plasma displays strong residual interactions through soft gluonic modes, which cannot be treated perturbatively, and which influence different quantities in different ways. In particular, any constituents themselves are objects dressed by non-perturbative interactions. This gives a consistent explanation to the various observed features: the equation of state, susceptibilities and fermionic correlators are dominated by hard modes, but significant deviations from ideal gas behavior are still present. Gluonic correlators, on the other hand, are dominated by soft modes and entirely off their leading perturbative predictions. An ideal gas is established only at asymptotically high temperatures.
The study of real time dynamics in the plasma is a much more difficult problem, for which at present no formal approach exists within lattice gauge theory. A step towards inspecting such physics could be presented by means of maximum entropy methods. However, at this stage such calculations have still a qualitative character, in trying to establish a working methodology.
A large section of this paper was devoted to another new development, allowing to study QCD also at small baryon density, as it is explored by heavy ion experiments at RHIC and LHC. Again, the presented studies are still at an exploratory level, subject to the already mentioned systematic effects. However, rather different methods have led to a consistent dependence of the critical temperature on the chemical potential, giving confidence in their validity up to µ B ∼ 500 MeV. In this range the pseudo-critical line is found to be rather flat. There also is a prediction for a critical point to be cross checked in simulations to come.
Advances in lattice QCD have always been brought about by a mixture of new calculational schemes, new algorithms and refined analysis techniques. They have, however, also been related to progress in computational resources. With a new generation of machines in the Teraflops-range going into operation soon, one has all reasons to hope that part of the mentioned systematic deficiencies will be reduced and simulations at parameter values much closer to the physical ones will be realized. We hope to have also shown that, by combining analytic and numerical techniques, insight can be gained and new paths may show up, which might be needed to venture fully into the real time regime.