QCD thermodynamics on the lattice

A remarkable progress has been made in the understanding of the hot and dense QCD matter using lattice gauge theory. The issues which are very well understood as well those which require both conceptual as well as algorithmic advances are highlighted. The recent lattice results on QCD thermodynamics which are important in the context of the heavy ion experiments are reviewed. Instances of greater synergy between the lattice theory and the experiments in the recent years are discussed where lattice results could be directly used as benchmarks for experiments and results from the experiments would be a crucial input for lattice computations.


I. INTRODUCTION
A large part of the visible matter in our universe is made up of protons and neutrons, collectively called hadrons.
Hadrons are made up of more fundamental particles called quarks and gluons. The quantum theory for these particles is Quantum Chromodynamics(QCD). QCD is a strongly interacting theory and the strength of interaction becomes vanishingly small only at asymptotically high energies. Due to this reason, the quarks and gluons are not visible directly in our world and remain confined within the hadrons. Lattice gauge theory has emerged as the most successful non-perturbative tool to study QCD, with very precise lattice results available for hadron masses and decay constants which are in excellent agreement with the experimental values [1].
It is expected that at high enough temperatures that existed in the early universe, the hadrons would melt into a quark gluon plasma(QGP) phase. Signatures of such a phase have been seen during the last decade in the Heavy ion collision experiments at the Relativistic Heavy Ion Collider(RHIC), in Brookhaven National Laboratory This is particularly exciting for the lattice theory community which has been predicting such a phase transition since a long time [2]. The formation of the QGP phase occurs at temperatures near Λ QCD , where QCD is strongly interacting, which means lattice is the most reliable tool to understand the properties of the hot QCD medium. Over the past three decades, the lattice community has contributed significantly to the understanding of the physics of heavy ion experiments and strongly interacting matter under extreme conditions, in general. Lattice computations are entering into the precision regime, where lattice data can be directly used for interpreting the experimental results and set benchmarks for the heavy ion experiments at RHIC and at the ALICE facility in CERN. It is now generally believed that the hot and dense matter created due to the collision of two heavy nuclei at RHIC and ALICE, equilibrates within 1 fm/c of the initial impact. The equilibrated QGP medium then expands and cools down in the process, ultimately forming hadrons at the chemical freezeout. The evolution of the fireball from its equilibration till the chemical freezeout is described by relativistic hydrodynamics [3]. The QCD Equation of State(EoS) is an input for the hydrodynamic equations and lattice can provide a non-perturbative estimate of this quantity from first principles. The lattice data for the speed of sound in the QCD medium is also an important input for the hydrodynamic study, once bulk viscosity is considered.
In this article, I have selected the most recent results form lattice QCD thermodynamics that are relevant for the heavy ion phenomenology. I have tried to review the necessary background, but not attempted to provide a comprehensive account of the development of the subject throughout these years. I have divided this article into two major sections. The first section deals with QCD at finite temperature and zero baryon density, where lattice methods are very robust. I have given a basic introduction to the lattice techniques, and how the continuum limit is taken, which is essential to relate the lattice data with the real world experiments. I have discussed the current understanding we have of the nature of QCD phase transition as a function of quark masses, inferred from lattice studies. Subsequently the different aspects of the hot QCD medium for physical quark masses are discussed; the EoS, the nature and the temperature of transition and the behaviour of various thermodynamic observables in the different phases. In the study of thermodynamics, the contribution of the lighter u, d and s quarks are usually considered. The effect of heavier charm quarks on QCD thermodynamics is discussed in this section, in view of their relevance for the heavy ion experiments at LHC, where hydrodynamic evolution is expected to set in already at temperatures about 500 MeV and also for the physics of early universe. The relevance of chiral symmetry for the QCD phase diagram and the effects of chiral anomaly are discussed in detail. The chiral anomaly is believed to have an important role in shaping the phase diagram and several lattice studies in the recent years are trying to understand its effect. It is a difficult problem and I have tried to compile the recent results and review the general understanding within the community, about how to improve upon them.
The second section is about lattice QCD at finite density, where there is an inherent short-coming of the lattice algorithms due to the so-called sign problem. A brief overview of the different methods used and those being developed by the lattice practitioners to circumvent this problem, is given. It is an active field of research, with a lot of understanding of the origin and the severity of this problem gained in recent years, which is motivating the search for its possible cure. In the regime where the density of baryons is not too large, which is being probed by the experiments at RHIC, lattice techniques have been used successfully to produce some interesting results. One such important proposal in the recent time, is the first principles determination of the chemical freezeout curve using experimental data on the electric charge fluctuations. This and the lattice results on the fluctuations of different quantum numbers in the hot medium and the EoS at finite baryon density are discussed in detail. An important feature of the QCD phase diagram is the possible presence of a critical end-point for the chiral first order transition. Since critical end-point search is one of the main objectives at RHIC, I have reviewed the current lattice results on this topic. The presence of the critical end-point is still not conclusively proven from lattice studies. It is a very challenging problem and I mention about the further work in progress to address this problem effectively. Fermions with exact chiral symmetry on the lattice are important in this context. I have discussed the recent successful development to construct fermion operators that have exact chiral symmetry even at finite density which would be important for future studies on the critical end-point. The signatures of the critical end-point could be detected in the experiments if the critical region is not separated from the freezeout curve. It is thus very crucial to estimate the curvature of the critical line from first principles and I devote an entire subsection to discuss the lattice results on this topic.
I apologize for my inability to include all the pioneering works that have firmly established this subject and also to review the extensive set of interesting contemporary works. For a comprehensive review of the current activity in lattice thermodynamics, at finite temperature and density, I refer to the excellent review talks of the Lattice conference, 2012 [4,5].

II. QCD AT FINITE TEMPERATURE ON THE LATTICE
The starting point of any thermodynamic study is the partition function. The QCD partition function for N f flavours of quarks in the canonical ensemble is given as, (1) where D f is the fermion operator for each flavour of quark f . U µ is the gauge link defined as U µ (x) = exp(−ig x A µ (x ′ )dx ′ ) in terms of gauge fields A µ , which are adjoint representation of the SU (3) color group and g is the strength of the gauge coupling. S G is the gluon action in Euclidean space of finite temporal extent of size denoted by the inverse of the temperature of the system, T . Lattice QCD involves discretizing the spacetime into a lattice with a spacing denoted by a. The volume of the lattice is given as V = N 3 a 3 , where N are the number of lattice sites along the spatial directions and the temperature being T = 1/(N τ a), where N τ are the number of sites along the temporal direction. The lattice is usually denoted as N 3 × N τ . The gluon action and the fermion determinant are discretized on the lattice. The simplest gluon action, known as Wilson plaquette action is of the form, where U µ,ν (x) is called a plaquette. The naive discretization of the continuum Dirac equation on the lattice results in the fermion operator of the form, where in each of the expressions the site index runs from x = 1 to N 3 × N τ . The discretization of the gluon and fermion operators are not unique and there are several choices which give the correct continuum limit. Usually discretized operators with small finite a corrections are preferred. Reducing a-dependent corrections by adding suitable "irrelevant" terms in the Renormalization Group(RG) sense, is known as improvement of the operator. Another issue related to the discretization of the fermion operator is called as the "fermion doubling problem". It arises because the naive discretization of the continuum fermion operator introduces extra unphysical fermion species called the doublers. The existence of the doublers can be traced back to a No-Go theorem [6] on the lattice which states that fermion actions which are ultra-local, have exact chiral symmetry, and have the correct continuum limit, cannot be free from the doublers. Doublers are problematic since in the continuum limit we would get a theory with 16 fermion species and QCD with 16 flavours is very close to the upper bound of the number of flavours beyond which the asymptotic freedom is lost. It is thus important to ensure that the discrete fermion operator should be free of the doublers. In order to do so the chiral symmetry is explicitly broken on the lattice, like for the case of Wilson fermions [7], or only a remnant of it is preserved for the staggered fermions [8]. The staggered fermion discretization retains the doubling problem in a milder form. In the continuum limit, the staggered fermion determinant would give contribution of four degenerate fermion species or tastes. However on a finite lattice, there is a considerable mixing among the tastes so a simple fourth root of the determinant would not yield the contribution of a single fermion flavour. This is called the rooting problem. The severity of rooting problem can be minimized by choosing either the stout-smeared staggered quarks [9] or the Highly Improved Staggered Quarks(HISQ) [10]. Other improved versions of staggered fermions used for QCD thermodynamics are the p4 and asqtad fermions [11,12]. Only the overlap [13] and the domain wall fermions [14] have exact chiral symmetry on the lattice at the expense of breaking the ultra-locality condition of the Nielsen-Ninomiya No-go theorem. As a result overlap and domain wall fermions are much more expensive to simulate compared to the staggered and the Wilson fermions. For QCD thermodynamics, the staggered and to some extent the Wilson fermions are favourites, with very high precision data available with improved versions of staggered quarks [15,16]. With the advent of faster computing resources and smarter algorithms, even large scale simulations with chiral fermions are becoming a reality [17][18][19][20].
With the choice of a suitable gauge and the fermion operators on the lattice, different physical observables are measured on statistically independent configurations generated using suitable Monte-Carlo algorithms. To make connection with the continuum physics, one needs to take the a → 0 limit of the observables measured on the lattice. The gauge coupling is related to the lattice spacing through the beta-function and the continuum limit, in turn, implies g → 0. In the space of coupling constants and the fermion masses, the continuum limit is a second order fixed point and the approach to the fixed point should be done along the correct RG trajectory or the lines of constant physics. The line of constant physics is defined by setting the mass of hadrons on the lattice to the continuum values, at each value of the coupling constant. The number of such relations required depends on the number of fermion flavours. To relate the lattice hadron masses to their experimental values, one has to define a scale to express the lattice spacing a, in terms of some physical units. There are two often used methods in QCD to set the scale, using the quantities r 1 and the kaon decay constant f K . The r 1 scale is defined from the quark-antiquark potential Vq q (r), as, On the lattice one measures Vq q (r) and r 1 is extracted from it using a suitable fit ansatz for the potential. To quantify the value of r 1 in physical units, one uses either the pion decay constant or the splitting of energy levels of bottom mesons to set the lattice spacing [21]. Advantages of this scale is that it is not sensitive to fermion discretization effects and to the choice of quark masses that defines the line of constant physics. However, the accurate determination of the potential requires very good statistics. One can also set the scale by choosing the f K measured on the lattice to its physical value. The f K is known with very high accuracy from the experiments. Once the line of constant physics is set, one has to take care of the finite size and lattice spacing effects such that the continuum extrapolation is correctly performed. To minimize such corrections, the correlation length which is given by the inverse of the mass of the lowest excitation of the system should be much larger that the lattice spacing but sufficiently smaller than the spatial size. Also for thermodynamics, it is crucial to minimize finite volume corrections which is ensured for the choice ζ ≥ 3, where ζ = N/N τ .
To characterize different phases one needs to define a suitable order parameter which depends on the symmetries of the theory. In the limit of infinitely heavy quark masses QCD is just a pure gauge theory with an exact order parameter, the expectation value of the Polyakov loop given as, The phase transition from a phase of confined colour degrees of freedom to the deconfined regime of free gluons, is of first order and is established very firmly from lattice studies [22]. The corresponding transition temperature is T c (pure-gauge)= 276(2) MeV [23] using the string tension, √ σ, value to be 425 MeV, to set the scale . If the quarks are massless, the QCD partition function with N f quark flavours has an exact SU (N f ) ⊗ SU (N f ) chiral symmetry. At some temperature, there is a phase transition from a chiral symmetry broken phase to the symmetry restored phase, characterized by the order parameter called the chiral condensate, The phase transition in the chiral limit for N f = 3, is expected to be of first order and there are several lattice results supporting this [24]. For N f = 2, the lattice results are contradictory with some claiming a first order transition [25] whereas recent results showing that the second order transition is also a possibility [26]. The current status of N f = 2 QCD phase transition in the chiral limit would be discussed again in a later subsection. For any finite value of quark masses, however there is no unique order parameter and no sharp phase transition is expected but only a gradual crossover. Based on effective field theories with same symmetries as QCD, using universality arguments and renormalization group inspired techniques, a schematic diagram of different phases of QCD as a function of quark mass is summarized in the famous "Columbia plot" [27]. The first order regions in the quenched and the chiral limits are separated from the cross-over region by second order lines belonging to the Z(2) universality class. These boundaries are schematic, though, and it is important to estimate the precise location of the physical point in this diagram. Lattice studies over the years have helped to redraw the boundaries more quantitatively. A latest version of the "Columbia plot" is shown in figure 1. With the high precision lattice data with physical light and strange quark masses, it is now known that the QCD transition in our world is a crossover [28][29][30]. The boundary of the first order region in the upper right hand corner of figure 1 is fairly well known [31]. The extent of the first order region in the bottom left hand is now believed to be small and much far away from the physical point [32,33]. However the extent of the Z(2) line in the left hand corner is still not well established; it can either continue along the m u,d = 0 axis to the m s → ∞ corner or end at a tricritical point. A better understanding of this issue is currently underway. The key to the resolution of this issue is to understand the effects of chiral anomaly through rigorous lattice computations. Since the light u,d-quark masses are much smaller than Λ QCD , the QCD action has an approximate SU (2) × SU (2) × U B (1) symmetry with an additional classical U A (1) symmetry broken explicitly by quantum effects. This is known as the U A (1) anomaly [34,35]. At zero temperature, the magnitude of this anomaly is related to the instanton-density. If the magnitude of this anomaly is temperature independent, the phase transition along the m u,d = 0 axes has to be of second order, belonging to the O(4) universality class [36]. This would mean that the Z(2) line has to end at a tri-critical point characterized by the strange quark mass, m tric s . The difference between the physical and tri-critical mass for the strange quark is not yet known with a good precision. In the following subsections, the lattice results for the QCD EoS for physical quark masses are discussed, which is an input for the hydrodynamics of the QGP medium. The current results on the pseudo-critical temperature, the entropy density, the speed of sound are also shown. All the results are for 2 + 1 flavour QCD, i.e., two light degenerate u and d quarks and a heavier strange quark mass. The effect of the heavy charm quarks on the thermodynamic quantities is also highlighted. At the end of this section, I touch upon the N f = 2 QCD near the chiral limit and the effects of the U A (1) anomaly for QCD thermodynamics.  (ε-3p)/T 4 HISQ/tree:

FIG. 2:
The results for the trace anomaly using the HISQ action for low (left panel) and high(right panel) temperatures for lattice sizes with temporal extent Nτ and spatial size 4Nτ , from [39]. Also in the right panel, the HISQ results are compared to the results using p4 fermions, which has an improved behaviour at high temperatures and to the continuum perturbation theory results at 1-loop(solid line) and 2-loop(dashed line)for the trace anomaly. The stout data are the continuum estimates from the Nτ = 6, 8, 10 data in [16].

A. Equation of state
The Equation of State(EoS) is the relation between the pressure and energy density of a system in thermal equilibrium. For estimating the QCD EoS, the most frequently used method by the lattice practitioners is the integral method [37]. In this method, one first computes the trace anomaly I(T ), which is the trace of the energy-momentum tensor. This is equal to the quantity ǫ − 3p where ǫ is the energy density of the system and p is the pressure. Moreover it is related to the pressure of the system through the following relation So if I(T ) is known, the pressure can be computed by integrating I(T ) over a range of temperature, with the lower value of temperature chosen such that the corresponding value of pressure is vanishingly small. The trace anomaly is related to the chiral condensate and the gluon action as, where the subscript zero denotes the vacuum expectation values of the corresponding quantities. The subtraction is necessary to remove the zero temperature ultraviolet divergences and the vacuum expectation values are usually computed on a lattice with number of sites (N τ ) 0 in the temporal direction, equal to the corresponding spatial number of sites, N . The subtraction is an unavoidable expense of this method. A new idea of deriving thermodynamic observables from cumulants of momentum distribution has emerged, where the vacuum subtraction is not required [38] and it would be interesting to check the application of this method in QCD. Also one needs to know the functional dependence of the inverse of QCD coupling constant β and the quark masses with the lattice spacing a along the line of constant physics. On the lattice, I(T ) is known only for a finite number of temperature values. The pressure computed by the numerical integration of the I(T ) data, has errors both due to statistical fluctuations and systematic uncertainties involved in the numerical interpolation of the data. The results for the trace anomaly are available for different lattice discretizations of the fermions. For staggered quarks there are two sets of results, one from the HotQCD collaboration using HISQ discretization [39,40] and the other from the Budapest-Wuppertal collaboration using stout smeared staggered quarks [16,41]. These results are compiled in Figures 2, 3. For the HISQ results, the bare lattice parameters are fixed by setting the lowest strange pseudo-scalar meson mass to its physical value at about 686 MeV and m π = 160 MeV, which defines the line of constant physics. The kaon decay constant f K = 156.1 MeV or alternatively the r 1 = 0.3106 fm from the static quark potential is used to set the scale. The corresponding parameters for the stout smeared quarks are m π = 135 MeV, m K = 498 MeV and the kaon decay constant. From figure 2, it is evident that there is a good agreement between the two sets of results for T < 180MeV and also for high enough temperatures T > 350 MeV. The stout continuum results in the figure were obtained extrapolation with the N τ = 6, 8, 10 data from Ref. [16]. In the intermediate temperature range there is some discrepancy, specially the peaks of the interaction measure do not coincide for these two different discretization schemes, which may be due finite lattice spacing effects. However the HISQ N τ = 12 data is inching closer to the stout results in this regime. The recent continuum stout results, obtained from continuum extrapolation of the new N τ = 12 data in addition to the older data are consistent with the HISQ results, with the peak position shifting to 200 MeV(left panel of figure 3). There is also a good agreement of the HISQ and stout data with the trace anomaly obtained from the Hadron Resonance Gas(HRG) model for T < 140 MeV and with the resummed perturbation theory results at high temperatures. Using the N τ = 6, 8 data which is available upto temperatures of 1000 MeV, a continuum extrapolation of the stout data was performed, the result of which is shown in right panel of figure 3. For this entire range of temperature, there is a useful parameterization characterizing the trace anomaly [16] with the following parametric form, where the best fit parameters are This parametric form could be a useful input for the hydrodynamical simulations, which usually uses the lattice EoS before hadronization and that from the HRG after the freezeout of hadrons. There are lattice results for the EoS using alternative fermion discretizations, the Wilson fermions. The WHOT-QCD collaboration has results for 2 + 1 flavours of improved Wilson fermions [42] with the physical value of strange quark mass but a large pion mass equal to 0.63m ρ . The tmfT collaboration has results for 2 flavours of maximally twisted Wilson fermions [43] with m π > 400MeV. Both these results are compiled in figure 4. These are in rough qualitative agreement with the staggered fermion data, specially the peak for the WHOT-QCD data occurring at 200 MeV is consistent with the HISQ and stout results. A more quantitative agreement at this stage is difficult, since the pion masses for the Wilson fermions are much larger than the physical value.

B. The pseudo-critical temperature
We recall that the QCD transition, from a phase of color singlet states to a phase of colored quantum states is an analytic crossover, for physical quark masses. This is fairly well established by now from lattice studies using two different approaches. One approach is to monitor the behaviour of the thermodynamic observables in the transition region for physical values of quark masses while the other is to map out the chiral critical line as a function of light quark mass [44]. The absence of a sharp phase transition implies that there is no unique transition temperature but only different pseudo-critical temperatures corresponding to different observables. There is no order parameter but the observables like the renormalized Polyakov loop, L R has a point of inflexion across the crossover region. Another observable relevant in the crossover regime is the renormalized chiral condensate, which has been defined [45] in the following manner to take into account the multiplicative renormalization as well additive ones due to a finite bare  quark mass, The normalized chiral susceptibility χ R for the light quarks defined as, is a good observable as well. Both L R and ∆ l,s (T ) have a point of inflexion at the pseudo-critical temperature and χ R has a smooth peak. From the continuum extrapolated data of the stout-smeared staggered fermions, the pseudo-critical temperatures corresponding to these observables for physical quark masses are, The data for L R and ∆ l,s with the HISQ disretization, is shown in figure 5. These are for lattice of size N τ × (4N τ ) 3 . The HISQ data are in good agreement with the continuum extrapolated stout-smeared staggered results from [46]. The fact that the rise of L R is more gradual than the corresponding rise of ∆ l,s signals that the crossover is more likely influenced by the chiral symmetry restoration. Previous scaling studies of the renormalized chiral condensate with the p4-staggered quarks showed that the physical light quarks already approximate the O(4) critical behaviour of the chiral quarks [26]. Using the O(4) scaling of the renormalized chiral condensate, the T c obtained for HISQ quarks through chiral and continuum extrapolation is 154 ± 9 MeV. This value is in excellent agreement with the stout result, implying that the continuum extrapolation done with the staggered fermions is quite robust.

C. Comparing results for different fermion discretizations
The results for the EoS and the pseudo-critical temperature discussed so far, have been obtained using different improved versions of the staggered quarks. For these fermion species, the so called "rooting" problem may alter the continuum limit due to breaking of the U A (1) anomaly [47] though some other work refutes this claim [48]. It is important to check the effects of the rooting procedure on the continuum extrapolation of finite temperature observables. The Budapest-Wuppertal collaboration has recently compared the continuum extrapolated results for different observables using the Wilson and staggered fermions [49] as the former discretization does not suffer from the rooting problem. The scale for the Wilson fermions was determined using m Ω = 1672 MeV and the line of constant physics was set using m π /m Ω ∼ 0.3 and m K /m Ω ∼ 0.36. For the staggered quarks, the line of constant physics was set such that the ratios m π /m Ω and m K /m Ω are within 3% of the corresponding values for the Wilson fermions. This means that the pions are quite heavy with m π ∼ 540 MeV for both these discretizations. The continuum extrapolated results for L R and the renormalized chiral condensate are shown in figure 6. The continuum results for both these quantities are in good agreement for the whole range of temperature, implying that these two different fermion discretizations indeed have the correct continuum limit. In all these computations an improved ∆ l,s L ren (T) FIG. 5: The results for the subtracted chiral condensate(left panel) and the renormalized Polyakov loop(right panel) from the HotQCD collaboration, from [40]. These data are compared with the continuum results using stout smeared fermions, from [46].
Wilson operator was used, in which the dominant O(a) correction terms due to explicit breaking of chiral symmetry by these fermions were cancelled. It ensured that in both the studies the approach to the continuum limit was chosen to be the same. However at this large value of quark masses, the rooting problem may be mild enough to show any adverse effects and it would be desirable to perform a similar comparison at physical value of the quark masses. Since the effects of chiral symmetry persist in the cross-over region, it is important to compare the existing results for T c with those using fermions with exact chiral symmetry on the lattice. For the Wilson and the staggered action, even for massless quarks, the full SU (2)⊗SU (2) chiral symmetry is realized only in the continuum limit. For chiral fermions on the lattice, like the overlap or the domain wall fermions, the chiral and the continuum limits are disentangled, allowing us to understand the remnant effects of chiral symmetry in the cross-over region even on a finite lattice. However, lattice QCD with overlap fermions is computationally prohibitive [50] and currently better algorithms are being developed to simulate them with comparatively lesser effort [51]. The domain wall fermions have exact chiral symmetry only when the extent of the fifth dimension, N 5 , of the five dimensional lattice on which these fermions are defined, is infinite. For smooth gauge fields, the chiral symmetry violation on a finite lattice is suppressed as an exponential of N 5 but the suppression could be much slower, as 1/N 5 for rough gauge configurations in the crossover region. Better algorithms have been employed to ensure exponential suppression even for rough gauge fields [52]. The most recent results for the overlap fermions from the Budapest-Wuppertal collaboration [20] and the domain wall fermions from the HotQCD collaboration [52], are shown in figure 7. The renormalized chiral condensate for the overlap fermions are qualitatively consistent with the continuum staggered fermion results, even for small volumes and large pion masses of about 350 MeV around the crossover region. The lattice cut-off effects seem to be quite small for N τ = 8. The renormalized chiral condensate and the ∆ l,s for the domain wall fermions are shown in figure 7. The lattice size is 16 3 × 8 with the number of lattice sites along the fifth dimension taken to be 32 for T > 160 MeV and 48 otherwise and the pion mass is about 200 MeV. The lattice volume is comparatively small therefore these results do not show a sharp rise in the crossover region. With larger volumes, the rise in these thermodynamic quantities is expected to be much steeper. The value of T c estimated from the peak of the chiral susceptibility i.e the derivative of the chiral condensate is between 160-170 MeV which is consistent with the continuum results from the HISQ fermions. The renormalized chiral condensate for the overlap quarks is compared to the continuum extrapolated results using the stout smeared staggered quarks in the left panel, from [20]. In the right panel, the behaviour of different chiral condensates defined using the domain wall fermions are shown in the critical region, from [52].

D. The thermodynamical observables
Thermodynamic observables characterize the different phases across a phase transition. From the behaviour of these observables, one can infer about the degrees of freedom of the different phases and the nature of the interactions among the constituents. It was already known from an important lattice study that the pressure in high temperature phase of QCD showed a strong dependence on the number of quark flavours [53], signaling deconfinement of the quark and gluon degrees of freedom. Recent results for the pressure, entropy density and the speed of sound for QCD, using the stout-smeared staggered quarks are compiled in figure 8. Though in our world there is no real phase transition, the entropy density increases rapidly with temperature, again signaling the liberation of a large number of colour degrees of freedom. The entropy density for QCD is almost 20% off from the value of a free gas of quarks and gluons, even at temperatures about 1000 MeV. The deviation of the pressure of QGP computed at similar temperatures, from its free theory value is even more, close to about 25% of its value. Another observable that characterizes the different phases is the speed of sound c s . If QGP at high temperatures was qualitatively close to a strongly interacting conformal theory, then the speed of sound would be exactly 1/ √ 3. However the deviation from conformality is quite significant even at temperatures about T = 500 MeV which hints that the AdS-CFT inspired study of the QGP medium should be done with more care. The values of entropy density computed with different discretizations of FIG. 8: The entropy density, pressure and the speed of sound for the stout-smeared fermions as a function of temperature, from [16]. staggered fermions like the asqtad or the p4 fermions [54] show about 10% deviation from the free theory value at very high temperatures. The departure from AdS-CFT values is even more severe using these fermions. The stout results are about 10% lower than the corresponding asqtad and p4 results. This deviation is attributed to the fact that the latter discretizations have smaller cut-off effects at higher temperatures and would be more closer to the continuum results. The stout continuum values shown in the figure were obtained by averaging the N τ = 8, 10 data. A proper continuum extrapolation of the results for both the fermion discretizations is necessary for resolving the difference and for use of these values for the real world calculations. However, the lattice results with at least 10% off from the free theory values even at very high temperatures, implies that the QGP phase is strongly interacting, more like a liquid rather than a gas of quarks and gluons, confirming the similar prediction from the RHIC experiments. For T < T c , the results for all these observables are in agreement with Hadron resonance gas model predictions.

E. Effects of charm quarks on the EoS
The effects of charm quarks to the pressure in the QGP phase was estimated sometime ago, using next-to leading order perturbation theory [55]. It was observed that the contribution of charm quarks become significant for temperatures T > 2T c . Preliminary data from the LHC already indicates that the charm quarks would thermalize quickly as the lighter quarks. It would then affect the EoS and thus the hydrodynamical evolution of the fireball formed at LHC energies. Lattice studies are important to quantify the contribution of charm to the EoS in the QGP phase. The first lattice studies were done by the RBC [56] as well as the MILC collaboration [57] with quenched charm quarks, i.e., by neglecting quantum fluctuations due to the charm quarks. The quenched charm results for the EoS differ from the 2+1 flavour results, already at 1.2 T c . Recent results from the Budapest-Wuppertal collaboration with dynamical charm quarks [41], however, show that the effects of charm quarks show up only around 300 MeV, more in agreement with the perturbative estimates. Both the approaches highlight the fact that the effects of charm quark should be considered for the EoS used as an input for the hydrodynamical evolution of the fireball at LHC energies, which may set in at T ∼ 500 MeV. It would be also important for the EoS of the Standard Model, important for the cosmological evolution in the early universe [58,59].  [57]. The lattice size is 24 3 × 6. In the right panel, the effects of dynamical charm quarks to the pressure is shown as a function of temperature, from [41]. The chiral phase transition for N f = 2 QCD is still not well understood from lattice studies, as was emphasized at the beginning of this section. Though the lattice results for 2 + 1 flavours with different fermion discretizations are in good agreement, the corresponding ones for the two light flavour case are still inconclusive. Two major approaches have been undertaken in the recent years to understand the order of this transition. One of them is to check the scaling properties of the order parameter. If the phase transition is indeed a second order one, then the order parameter would show O(4) scaling in the transition region. The second approach is to understand the effects of the U A (1) anomaly near in the phase transition. If the quantum fluctuations responsible for this U A (1) anomaly decrease significantly with temperature, it would result in the degeneracy of the masses of mesons of certain quantum numbers and a characteristic behaviour of the density of low lying eigenmodes of the fermion operator. I discuss the major lattice results using both these approaches, in the following paragraphs. Most of these approaches are hinting that the two flavour chiral phase transition may be a second order one.

Scaling analysis in the critical region
The order parameter that characterizes the chiral phase transition is the chiral condensate. A suitable dimensionless definition of the chiral condensate used in the lattice study by the BNL-Bielefeld collaboration [26] is, The additive ultraviolet divergences are not explicitly subtracted from the condensate and hence it is the bare value denoted by subscript b. This additive divergence would be included in the regular part and in the transition region, would be much smaller in magnitude than the singular part of M b . In the vicinity of the transition region, the order parameter can be written as, where f G is the universal scaling function, known from analysis of the O(N ) spin-models [60][61][62] with β and δ being the corresponding critical exponents. The quantities h and t are dimensionless parameters that determine the deviations from the critical point and are defined as, with T c,0 being the transition temperature in the chiral regime i.e, for h → 0 and h 0 and t 0 are non universal constants.
One of the choice of the regular part of the order parameter used in the lattice study is, where one assumes that the regular part is an analytic function of the relevant parameters around the transition point. The BNL-Bielefeld collaboration used an improved variety of the staggered quarks, called the p4 quarks, to compute the order parameter defined in Eq. (12) and χ m , its derivative with respect to m l for different values of the light quark masses, m l . The strange quark mass was fixed at its physical value. These quantities were fitted to the functional form given in Eq. (13) and its derivative respectively. The scaling analysis was done for a fixed lattice of size N 3 × 4, so the order parameter and its derivatives are expected to have an O(2) scaling in the chiral regime since the fermion discretization only retains a remnant of the continuum O(4) symmetry group. From the plots for the order parameter in the left panel of figure 10, it is evident that for m l /m s = 1/80 the phase transition is indeed a second order one with O(2) critical exponents, though O(4) scaling cannot be ruled out completely with the current precision available. In the scaling regime, the variable M b /h 1/δ should be a universal function of t/h 1/βδ . In the right panel of figure 10, the scaled chiral condensate is seen to be almost universal for m l /m s < 1/20, which provides a hint that even for the physical quark masses there is a remnant effect of the chiral symmetry. The crossover transition for 2+1 flavour QCD should be sensitive to the effects of chiral symmetry and therefore also to the effects of the U A (1) anomaly.

The effects of UA(1) anomaly
The QCD partition function breaks U A (1) symmetry explicitly. However its effects varies with temperature since we know that at asymptotically high temperatures, we approach the ideal Fermi gas limit where this symmetry is restored. It is important to understand the temperature dependence of U A (1) breaking near the chiral phase transition. If U A (1) breaking is significantly reduced from that at zero temperature, one would then claim that the symmetry is effectively restored. This would result in the degeneracy of the mass of the isospin triplet pseudoscalar(pion) and scalar (delta) mesons. The order parameter for such an effective restoration is the quantity defined as, and the order parameter for the restoration of the chiral symmetry is the chiral condensate. These quantities are also related to the fundamental theory through the density of eigenvalues, ρ(λ) of the Dirac operator as, Different scenarios that could lead to different functional behaviour of ρ(λ) were discussed in detail in Ref. [52]. I summarize the arguments below, • From dilute instanton gas approximation: ρ(λ, m) = c 0 m 2 δ(λ) ⇒ ψ ψ ∼ m and χ π − χ δ ∼ 2 • Analyticity of ρ(λ, m) as a function of λ and m when chiral symmetry is restored. To the leading order In fact to understand the effect of anomaly, it is desirable to use fermions with exact chiral symmetry on the lattice. The overlap and the domain wall fermions are such candidates, for which the chiral anomaly can be defined. Indeed, the overlap fermions satisfies an exact index theorem on the lattice [63]. A recent study of the eigenvalue spectrum with the domain wall fermions from the HotQCD collaboration [64] seems to favour ρ(λ, m) = c 0 m 2 δ(λ) + c 1 λ, for the density of eigenvalues. This would imply that in the chiral limit, the U A (1) anomaly would still survive when the chiral symmetry is restored. This is also consistent with the behaviour of χ π − χ δ as a function of temperature, shown in the left panel of figure 11. At crossover temperature around 160 MeV, the χ π − χ δ is far from zero, implying that the effects of the anomaly may be large in the crossover region. A recent theoretical study [65] with the overlap fermions shows that in the chiral symmetry restored phase where ψ ψ = 0, the eigenvalue density in the chiral limit should behave as, which would imply that χ π − χ δ → 0 as m → 0. Moreover it is argued that if an operator is invariant under some symmetry transformation then its expectation value becoming zero would not necessarily imply that the symmetry is restored, whereas the converse is true [65]. This would mean that the observable χ π − χ δ may not be a good candidate to study the U A (1) restoration. Rather the equality of the correlators of the pion and delta meson could be a more robust observable to indicate the restoration of the U A (1) symmetry. Recent results from the JLQCD collaboration with 2 flavours of overlap fermions seems to indicate that the U A (1) may be restored near the chiral symmetry restoration temperature, making it a first order transition [66]. Two of their main results are compiled in figure 12. The correlators of the scalar mesons become degenerate at about 196 MeV and at the same temperature a gap opens up in the small eigenvalue region of the eigenvalue spectrum. T = 196 MeV is slightly above the transition temperature which is near about 177 MeV. For T = 177 MeV, there is no degeneracy between the scalar and the pseudoscalar correlators and the density of zero modes is finite implying that the chiral symmetry is broken, which means that the U A (1) changes rapidly near the phase transition. However the lattice size is 16 3 × 8, which is small enough to introduce significant finite volume and cut-off effects in the present results.
FIG. 12: In the left panel, the quark mass dependence of eigenvalue distribution for the overlap quarks are compared at different temperatures, from [66]. In the right panel, the degeneracy of the scalar and pseudoscalar mesons for overlap fermions are shown at a temperature of 192 MeV which is slightly higher than the corresponding pseudo-critical temperature, from [66].
With the chiral fermions, the fate of U A (1) in the crossover region is still undetermined and more work needs to be done for conclusive understanding of this issue. With Wilson and staggered quarks, the anomaly is recovered only in the continuum limit. For fine enough lattice spacings, one can however check the behaviour of the low lying eigenmodes and the meson masses for different quantum numbers, to understand the effects of the remnant U A (1) anomaly using these fermions. From the eigenvalue distribution of HISQ operator shown in the left panel of figure  13, it is evident that the effect of U A (1) still persists at T = 330 MeV [67]. The long tail in the low lying eigenmodes is not a finite volume artifact since it persists even for very large volumes. However, the data is quite noisy and more statistics is required for making a final conclusion. The screening masses for the mesons of different quantum numbers were obtained from lattice studies [68] with improved Wilson fermions(right panel of figure 13). In the transition region, the scalar and pseudoscalar mesons are not degenerate and an agreement seen only for temperatures above 1.2 T c . However the input quark masses are quite large compared to the physical values and more data is needed to take a final call. At present, the effects of quantum anomalies are not yet understood from lattice studies.

III. LATTICE QCD AT FINITE DENSITY
QCD with a finite number of baryons is relevant for the physics of neutron stars and supernovae. It is the theoretical setup for the heavy ion physics phenomena occurring at low center of mass energy, √ s, of the colliding nuclei. Some of these low √ s collisions are being investigated at the RHIC and to be probed further with the start of the heavy ion experiments at FAIR, GSI and NICA, Dubna. In fact, an interesting feature of the QCD phase diagram is the critical end-point related to chiral symmetry restoration. The existence of the critical point has important consequences on the QCD phase diagram and it is the aim of the extensive Beam Energy Scan(BES) program at the RHIC to search for it.
To explain these experimental results from first principles, we need to extend the lattice QCD formulation to include the information of finite baryon density. One of the methods is to work in a grand canonical ensemble. In such an ensemble, the partition function is given by where the chemical potential µ is the Lagrange multiplier corresponding to the conserved number density N , that commutes with the QCD Hamiltonian H QCD . N can be the baryon number or the net electric charge . The µ enters into the lattice fermion action as exp(±µa) factors multiplying the forward and backward temporal links respectively [69,70], referred to as the Hasenfratz-Karsch method. The naive fermion operator at finite µ on the lattice would be of the form, This is not a unique way of introducing µ and it could be also done in several different ways [71]. The lattice fermion determinant at finite µ, like in the continuum, is no longer positive definite since (21) and the interpretation of DU detD f (µ)e −SG as a probability weight in the standard Monte Carlo simulations is no longer well defined. This is known as the "sign problem". One may consider only the real part of the fermion determinant for Monte Carlo algorithms and generate configurations, by so called phase quenching. Once the partition function is known in the phase quenched limit, one can then use the reweighting techniques to generate the partition function of the full theory at different values of µ. The expectation value of the phase of the determinant, needed for re-weighting, at some finite µ is given as where ∆F is the difference between the free energy densities of the full and the phase quenched QCD. For two degenerate quark flavours, the phase quenched theory is equivalent to a theory with a finite isospin chemical potential [72] and ∆F is the difference of free energies of QCD with finite baryon(quark) chemical potential and that at an isospin chemical potential. These two theories are qualitatively quite different and the sign problem results in a very small overlap between these two theories. For isospin QCD, the charged pions are the lightest excitations and these can undergo a Bose-Einstein condensation for µ > m π /2. The difference between the respective free energies in this regime is quite large, leading to a severe sign problem. This is an algorithmic problem that can arise for any theory which has chiral symmetry breaking. A better understanding of the sign-problem has been achieved in the recent years with a knowledge of the regions in the phase diagram with severe sign problem and those where it is controllable [73,74]. There are several methods followed to circumvent this problem on the lattice, some of which are listed below, • Reweighting of the µ = 0 partition function [75][76][77] • Taylor series expansion [78,79] • Canonical ensemble method [80,81] • Imaginary chemical potential approach [82][83][84] • Complex Langevin algorithm [85][86][87] • Worm algorithms [88,89] • QCD on a Lefschetz thimble [90] The Taylor series method has been widely used in the lattice QCD studies in the recent years, which has lead to interesting results relevant for the experiments. One such proposal is the determination of the line of chemical freezeout for the hadrons in the phase diagram at small baryon density, from first principles lattice study. It was first proposed that cumulants of baryon number fluctuations could be used for determining the freezeout parameters [91] on the lattice. Last year, another interesting suggestion was made [92], where the experimental data on cumulants of electric charge fluctuations could be used as an input to compute the freezeout curve using lattice data. This and some other results are discussed in the subsequent subsections. Most of the results are obtained with improved versions of staggered fermions. It has been known that the rooting problem may be more severe at finite density [93]. It is thus important to explore other fermion formulations as well for lattice studies. Wilson fermions have been used but it is important to use chiral fermions, especially for the study of the critical point. I outline in the next subsection, the theoretical efforts in the recent years that have led to the development of fermion operators at finite density with exact chiral symmetry on the lattice which can be used for future lattice studies on the critical point.

A. Chiral fermions at finite density
The contribution of the U A (1) anomaly is believed to affect the order of the chiral phase transition at zero density and hence is crucial for the presence or absence of the critical point. If the anomaly is not represented correctly at finite density, it may affect the location of the critical point in the phase diagram, if it exists. Overlap fermions have exact chiral symmetry on the lattice, in the sense that the overlap action is invariant under suitable chiral transformations known as the Luscher transformations [94]. It can be further shown that the fermion measure in the path integral is not invariant under Luscher transformations, and its change gives the chiral anomaly. The index theorem, relating the anomaly to the difference between the fermion zero modes, can be proved for them [63]. Thus the overlap fermions have the properties analogous to the fermions in the continuum QCD. In the continuum, it is known that the anomaly is not affected in presence of a finite baryon chemical potential. It would be desirable to preserve this continuum property with the overlap fermions as well, such that the physical properties important for the existence of the critical point are faithfully presented on a finite lattice. Defining an overlap fermion action at finite chemical potential is non-trivial as the conserved currents have to be defined with care [95]. The first attempt to define an overlap fermion operator at finite density [96] was done in the last decade, and an index theorem at finite µ was also derived for them. However these overlap fermions did not have exact chiral symmetry symmetry on a finite lattice [97]. Moreover, the index theorem for them was µ-dependent, unlike in the continuum. Recently, overlap fermion at finite density has been defined from the first principles [98], which has exact chiral symmetry on the lattice [99] and preserves the µ-independent anomaly as well. A suitable domain wall fermion action has been also defined at finite density [99], which was shown to reproduce the overlap action in the appropriate limit. It would be important to check the application of these overlap and domain wall fermion operators at finite µ for future large scale QCD simulations.

B. Correlations and fluctuations on the lattice
The studies of fluctuations of the conserved charges are important to understand the nature of the degrees of freedom in a thermalized medium and the interactions among them [100,101]. The diagonal susceptibility of order n, defined as, measures the fluctuations of the conserved quantum number X. In a heavy-ion experiment the relevant conserved numbers are the baryon number B and electric charge Q. The strangeness S, is zero at the initial time of collision of heavy nuclei but strange quark excitations are produced at a later time in the QGP and is also believed to be a good quantum number. These fluctuations can be computed exactly on the lattice at µ = 0 from the quark number susceptibilities [102]. Continuum extrapolated results for the second order susceptibilities of baryon number, strangeness and electric charge exist for both HISQ [103]and stout smeared staggered quarks [104]. The fluctuations of baryon number are very well explained by the hadron resonance gas model for T < 160 MeV. However, the fluctuations of the strangeness are usually larger than the HRG values by about 20% in the freezeout region characterized by 160 ≤ T ≤ 170 MeV. The electric charge fluctuations, on the other hand, are smaller than the corresponding HRG values by 10% in the same region. The ratio of χ Q 2 /χ B 2 (µ = 0) ≃ 0.29 − 0.35 in the freezeout region. A first principle determination of this ratio is crucial, as it would allow us to relate the net baryon number fluctuations with the net proton number fluctuations, which is an observable in the heavy ion experiments [103]. At high temperatures, these fluctuations slowly approach the corresponding free theory value, with the continuum extrapolated data for the baryon number susceptibility showing about 20% deviation from the free theory value even at 2T c [103]. The data are in good agreement with resummed perturbation theory estimates at these temperatures [105,106] indicating that the QGP is still fairly strongly interacting even at temperatures around 2T c .
To relate to the results of the heavy ion experiments at a lower collision energy, √ s, one has to compute the fluctuations on the lattice at a finite value of µ. The most widely used lattice method to compute the susceptibilities at a finite value of quark chemical potential µ, is through the Taylor expansion of the corresponding quantity at µ = 0 ,e.g., The light and strange quark susceptibilities have been computed at finite but small densities from Taylor expansion, using asqtad staggered quarks [57] and the ratios of baryon number susceptibilities using the unimproved staggered fermions [91] in the region of interest for the RHIC experiments. All these ratios agree well with the estimates from the HRG model [91], the results for which are compiled in the right panel of figure 16. The ratios of susceptibilities serve as a good observable for comparing the lattice and the experimental data since these are free from the unknown quantities, like the volume of the fireball during freezeout [107]. The higher order susceptibilities χ n , for n > 4, are important even in the µ = 0 regime. In the chiral limit, it is expected that the fourth order baryon number susceptibility would have a cusp and the sixth order would diverge with O(4) scaling at the critical temperature. Even for physical quark masses, χ B 6 for QCD would show oscillations near the pseudo-critical temperature and χ B 8 would have negative values in the same region [108], quite contrary to the HRG predictions. Thus the signatures of critical behaviour could be understood by the careful study of these quantities already at µ ∼ 0, which is probed by the experiments at LHC [108].
Other important quantities of relevance are the off-diagonal susceptibilities. These defined as are a measure of the correlations between different quantum numbers and hence good observables to estimate the effects of interactions in the different phases of the QCD medium. It has been suggested that the quantity C BS = −3χ BS 11 /χ S 2 is a good observable to characterize the deconfinement in thermal QCD [109]. If the strangeness is carried by quark like excitations, the value of C BS would be identity and would be much smaller than unity in the phase where only the baryons and mesons carries the strangeness quantum number. Recent results from the HotQCD collaboration using HISQ action [103] show that C BS approaches unity very quickly at around 200 MeV implying that almost no strange hadrons survive in the QGP phase above T c . This is compiled in the left panel of figure 14. The HotQCD data is consistent with the corresponding continuum extrapolated data with the stout smeared fermions [104]. Also C BS is not sensitive to the sea strange quark masses for T > T c since the first partially quenched results [110] for this quantity are consistent with the full QCD results. The other important observable is the baryon-electric charge correlation.
In the confined phase, electric charge in the baryon sector is mainly carried by protons and anti-protons, therefore the correlation would rise exponentially with temperature if this phase could be described as a non-interacting gas consisting of these particles. At high temperatures, however, quark like excitations would be important and their masses being much smaller than the temperature, this correlation would fall to zero. From the behaviour of the continuum extrapolated HISQ data for χ BQ 11 , compiled in the right panel of figure 14, it is evident that near the pseudo-critical temperature there is a change in the fundamental properties of the degrees of freedom of the medium, with quark like excitations dominating at 1.5T c . FIG. 14: The HISQ data for CBS(left panel) and χ BQ 11 /T 2 (right panel) as a function of temperature, from [103].

C. The freezeout curve from lattice
To relate the results from heavy ion experiments with the lattice data, it is crucial to map the center of mass energy of the colliding nuclei in the heavy ion collisions, √ s, to the corresponding point in the T − µ B plane of the QCD phase diagram. This is called the freezeout curve. Phenomenologically, the freezeout curve is obtained from a particular parameterization of the HRG model obtained through fitting the experimental data on hadron abundances [111]. At chemical freezeout, the chemical composition of the baryons gets frozen, meaning that the inelastic collisions between these species become less probable under further cooling of the system. However the systematic uncertainties in determining the hadron yields are not taken into account in the phenomenological determination of the freezeout curve. Recent work by the BNL-Bielefeld collaboration shows how lattice techniques can provide first principle determination of the freezeout curve through suitable experimental observables [92] . As emphasized in the last subsection, the ratios of susceptibilities are believed to be good observables for comparing the lattice and the experimental data. Two such observables proposed in Ref. [92] are, where M X , σ X , S X denotes the mean, variance and the skewness in dimensionless units for the conserved quantum number X. These observables are chosen because these have are odd and even functions of µ B , allowing us to independently determine T and µ B from these two quantities. The quantum number X can either chosen to be the net electric charge Q or the net baryon number B. In the experiments one can only measure the proton number fluctuations and it is not clear whether the proton number fluctuations could be a proxy for the net baryon fluctuation [112]. It was thus suggested, that the ratios of net charge fluctuations would be a better observable to compare with the experiments. Once the R Q 31 is known from experiments, one can determine the freezeout temperature T f from it by comparing with the continuum extrapolated lattice data. Analogously, one can obtain the µ B at freezeout from comparison of the R Q 12 data. In the left panel figure 15, the results for R Q 31 are shown as a function of temperature. It is evident that the first order correction to the value of the ratio is within 10% of the leading order value for µ B /T < 1.3 and in the freezeout region i.e, T > 140 MeV. From the leading order results of R Q 31 one can estimate the freezeout temperature. For √ s in the range of 39-200GeV currently probed in the Beam Energy Scan(BES) experiment at RHIC, the freezeout temperature from the HRG parameterization of the hadron multiplicities is about 165 MeV. At this temperature, the ratio R Q 31 calculated from the HRG model is quite larger than the lattice estimate which would mean that the freezeout temperature estimated from lattice data would differ from the model results by atleast 5%. Similarly, if R Q 12 is known from the experiments, µ B can be accurately estimated and is expected to be different from the current HRG estimates. This is not very surprising because the freezeout of the fluctuations happens due to diffusive processes and is due to a different mechanism from the freezeout of hadrons due to decreasing probability of inelastic collisions. Another question that was addressed in this work was how relevant are the other parameters like µ S and µ Q for the phase diagram and the freezeout curve. It was seen that µ S and µ Q are significantly smaller than µ B and the ratios of these quantities have a very small µ B dependence in the entire temperature range of 140-170 MeV relevant for the freezeout studies. It signifies that the relevant axes for the phase diagram are indeed T and µ B and these two parameters are sufficient for characterizing the freezeout curve. shown in the yellow band, is compared to its NLO value denoted by the blue band in the continuum limit. In the right panel, the ratios of µQ and µS with respect to µB are compared with the HRG model predictions at different temperatures. Both figures are from Ref. [92].

D. Physics near the critical point
It is known from models with the same symmetries as QCD, that the chiral phase transition at T = 0 and finite µ is a first order one. At zero density and high enough temperatures, QCD undergoes a crossover from the hadron to the QGP phase. By continuity, it is expected that the first order line should end at a critical end-point in the phase diagram [113]. The determination of its existence from first principles lattice computation has been quite challenging and the currently available lattice results are summarized in the left panel of figure 16. These are all obtained using staggered fermions. The first lattice study on the critical point was done using reweighting technique. Configurations generated at the critical value of the gauge coupling for µ B = 0 were used to determine the partition function at different values of T and µ B using two-parameter reweighting [76]. By observing the finite volume behaviour of the Lee-Yang zeroes of the partition function, it was predicted that for 2+1 flavour QCD, there is a critical end-point at T E = 160(4) MeV and µ B = 725(35) MeV. In this study the light quark was four times heavier than its physical value. Reducing the light quark mass, shifted the critical end-point to µ B = 360(40) MeV with T E = 162(2) remaining the same [114]. However, this result is for a rather small lattice of size 16 3 × 4 and is expected to change in the continuum limit and with larger volumes. Reweighting becomes more expensive with increasing volume of the lattice, so going to a larger lattice seems difficult with this method.
The other results for the critical point were obtained using the Taylor series method. In this method, the baryon number susceptibility at finite density is expanded in powers of µ B /T as a Taylor series as shown in Eq. (24), for each value of temperature. The baryon number susceptibility is expected to diverge at the critical end-point [115], so the radius of convergence of the series would give the location of the critical end-point [79]. However on a finite lattice, there are no divergences but the different estimates of the radius of convergence given as should all be positive and equal within errors at the critical end-point. Currently, the state of the art on the lattice are estimates of baryon number susceptibilities upto χ B 8 . This gives five different independent estimates of the radius of convergence upto r 6 , which were shown to be consistent within errors for N τ = 4, 6, 8 at T E = 0.94(1)T c [116][117][118]. The radius of convergence after finite volume correction is µ B /T E = 1.7(1) [118], which means µ B = 246(15) MeV at the critical end-point if we choose T c = 154 MeV. The input pion mass for this computation is about 1.5 times the physical value and could affect the final coordinates of the end-point. Moreover the different estimates for the radius of convergence r n in Eq. (27), agrees with each other for asymptotically large values of n and one might need to check the consistency of the results with the radii of convergence estimates beyond r 6 . Hints of the critical end-point were also obtained [119] using a different fermion discretization and a different methodology as well. Working with the canonical ensemble of improved Wilson fermions, the presence of a critical point was reported at T E = 0.925(5)T c and µ B /T c = 2.60 (8). This is a very preliminary study though, with a small lattice volume and a very heavy pion mass of about 700 MeV.
Though there is growing evidence in support for the existence of the critical end-point, the systematics for all these lattice studies are still not under control. It would be desirable to follow a different strategy to determine its existence. The alternate method suggested [120] was to determine the curvature of the surface of second order chiral phase transitions as a function of the baryon chemical potential µ B . If the chiral critical surface bends towards larger values of m u,d with increasing baryon chemical potential and for a fixed value of the strange quark mass, it would pass through the physical point, ensuring the existence of a critical end-point. However if the curvature is of opposite sign, the chiral critical end-point would not exist. For lattice size of 8 3 × 4, the critical value of the light quarks was estimated upto O(µ 4 B ) [121], with the strange quark mass fixed at its physical value. The leading value of the curvature has the same sign even for a finer lattice of extent N τ = 6 [122]. These studies show that the region of first order transition shrinks for small values of µ B , which would mean that the critical point may not exist in this regime of µ B . However for larger values of µ B , the higher order terms could be important and may bend the chiral critical line towards the physical values of quark masses. The finite cut-off effects are still sizeable and it is currently premature to make any definite predictions in the continuum limit with this method.  [118]. The magenta solid circle, box and star denote the Nτ = 4, 6, 8 data respectively for 2 flavours of staggered quarks [116][117][118] and the open circles denote Nτ = 4 data for 2+1 flavours obtained with reweighting techniques [76,114]. In the right panel, the ratio of the third and the second order baryon number susceptibility is plotted as a function of √ s relevant for the RHIC and LHC experiments and compared with the HRG model data, from [91].
It is equally important to understand the possible experimental signatures of the critical point. The search of the critical end-point is one of the important aims for the extensive BES program at RHIC. In a heavy ion experiment, one measures the number of charged hadrons at the chemical freezeout and its cumulants. During the expansion of the fireball, the hot and dense QCD medium would pass through the critical region and cool down eventually forming hadrons. If the freezeout and the critical regions are far separated, the system would have no memory of the critical fluctuations and the baryon number susceptibility measured from the experiments could be consistent with the predictions from thermal HRG models which has no critical behaviour. If the freezeout region is within the critical region, the critical fluctuations would be larger than the thermal fluctuations. It is thus important to estimate the chiral critical line for QCD from first principles. The curvature of the chiral critical line has been estimated by the BNL-Bielefeld collaboration [123], by extending the scaling analysis of the dimensionless chiral condensate M b outlined in Section II F 1 for finite values of baryon chemical potential, using Taylor series expansion. The corresponding scaling variables at finite µ B are The quantity M b can be expanded as a Taylor series in µ B /3T as, where χ m,B is the mixed susceptibility defined as χ m,B = T 2 ms ∂ 2 M b /∂(µ B /3T ) 2 computed at µ B = 0. In the critical region, it would show a scaling behaviour of the form, The universality of the scaled χ m,B data is clearly visible in the right panel of figure 17, both for p4 staggered quarks on N τ = 4 lattice with mass ratios of light and strange quarks varying from 1/20 to 1/80 and with HISQ discretization on a 32 3 × 8 lattice with the mass ratio fixed at 1/20. The fit of the complete lattice data set to the scaling relation for χ m,B , gave the value of κ B = 0.00656 (66). At non-vanishing µ B , the phase transition point is located at t = 0, which implies that the critical temperature at finite density can be parameterized as This estimate of the curvature is about three times larger than the corresponding prediction from the hadron resonance gas model. It would be interesting to compare the curvature of the freezeout line computed on the lattice with that of the critical line, once the experimental data for the electric charge cumulants are available.
Another complimentary study about the fate of the critical region at finite density was done by the Budapest-Wuppertal group [124]. It was suggested that if the critical region shrinks with increasing µ B , it would imply that one slowly converges to the critical end-point. The width of the critical region was measured from two different observables, the renormalized chiral condensate and the strange quark number susceptibility. Stout smeared staggered quarks were employed and the continuum limit was taken with the N τ = 6, 8, 10 data. The results are summarized in the left panel of figure 17. From the plots, it seems that the width of the crossover region does not change from its µ B = 0 value significantly for µ B < 500 MeV, which implies either that the critical end-point does not exist at all or is present at a higher value of µ B . The corresponding curvature measured for the light quark chiral condensate is 0.0066 (20) which is consistent with the result from the BNL-Bielefeld collaboration. The results indicate that the chiral pseudo-critical line and the phenomenological freezeout curve would separate apart at larger values of µ B and would be further away at the critical end-point.
It was noted that the higher order fluctuations are more strongly dependent on the correlation length of the system [125] and would survive even if the chiral and freezeout lines are far apart. It has been proposed [108], that the signature of the critical point can be detected by monitoring the behaviour of the sixth and higher order fluctuations of the electric charge along the freezeout curve. In the left panel, the width of the pseudocritical region for chiral condensate is shown as a blue curve and that for strange quark susceptibility is shown as a red curve, from [124]. In the right panel, the scaling of the mixed susceptibility is shown for different light quark masses and at the physical value of strange quark mass, from [123].

E. The EoS at finite density
The EoS at finite density would be the important input for understanding the hydrodynamical evolution of the fireball formed at low values of the collisional energy, at the RHIC and the future experiments at FAIR and NICA. It is believed that there is no generation of entropy once the fireball thermalizes [126]. In that case, as pointed out in Ref. [127], it is important to determine the EoS along lines of constant entropy per net baryon number, S/n B to relate the lattice results with the experiments. The isentrope, determined by a fixed value of S/n B that characterizes the evolution of the fireball, is S/n B ≃ 300 for RHIC experiments, at √ s = 200 GeV. For the future experiments at FAIR, the isentropes would be labeled by S/n B = 30 nearly as same as the early SPS experiments at CERN where S/n B ∼ 45. For two flavour QCD with p4 staggered quarks and with pion mass heavier than its physical value, it was already observed that the ratio of pressure and energy density showed little variation as a function of S/n B . The pressure and the energy density at finite µ are usually computed on the lattice as a Taylor series about its value at zero baryon density, as, P (µ l , µ s ) The formula is valid for two degenerate light quark flavours and a heavier strange quark. The coefficients χ ij are the quark number susceptibilities at µ = 0 and are non-zero for i + j=even. The corresponding expression for the trace anomaly is given as, The χ ij s can be also obtained from the coefficients b ij by integrating the latter along the line of constant physics. For 2+1 flavours of improved asqtad staggered quarks with physical strange quark mass and m l = m s /10, the interaction measure was computed upto O(µ 6 ) for two different lattice spacings (the left panel of figure 18). The interaction measure did not change significantly from the earlier results with heavier quarks and showed very little sensitivity to the cut-off effects along the isentropes [57]. However, it was observed that the light and the strange quark number susceptibilities change significantly from the zero temperature values along the isentropes. No peaks were found in the quark number susceptibilities at isentropes S/n b = 300, which led to the conclusion that the critical point may not be observed at the RHIC [57]. The EoS and the thermodynamic quantities were computed for physical values of quark masses by the Budapest-Wuppertal collaboration [128]. In this study, they set the values of the light quark chemical potentials such that µ l = µ B /3 and the strange quark susceptibility is µ s = −2µ l χ us 11 /χ s 2 to mimic the experimental conditions where the net strangeness is zero. The pressure and the energy density was computed upto O(µ 2 ). The ingredients that went into the computations were a) the near continuum values of the interaction measure data from the N τ = 10 lattice and b) the spline interpolated values of χ s 2 , χ us 11 for the range 125 < T < 400 MeV obtained using the continuum extrapolated data for χ s 2 , χ us 11 . It was observed, as evident from the right panel of figure 18, that the finite density effects along the RHIC isentropes are negligible consistent with the earlier work. However for isentropes given by S/n B = 30, the finite density effects become more important. The effect of truncation at O(µ 2 ) was also estimated on a reasonably large N τ = 8 lattice. It was observed, implying that the fourth and higher order terms need to be determined for even modest values of µ B in the Taylor series method. An independent study about the truncation effects of the Taylor series was performed in Ref. [129]. The derivatives of pressure were computed for two flavour QCD with staggered quarks at imaginary chemical potential. These derivatives are related to the successive terms of the Taylor coefficients of pressure evaluated at µ = 0. By fitting the imaginary µ data with a polynomial ansatz, these Taylor coefficients were obtained and compared with the exact values. It was observed that for T c ≤ T ≤ 1.04T c , at least the 8th order Taylor coefficient is necessary for a good fit. This highlights the necessity to evaluate higher order susceptibilities, beyond the currently measured eighth order in the studies of EoS or the critical end-point. New ideas to extend the Taylor series to higher order susceptibilities are evolving [129,130] and these should be explored in full QCD simulations.

IV. SUMMARY
As emphasized in the introduction, I have tried to compile together some of the important instances to show that the lattice results have already entered into the precision regime with different fermion discretizations giving consistent continuum results for the pseudo-critical temperature and fluctuations of different quantum numbers. The continuum result for the EoS would be available in very near future, with consistency already observed for different discretizations. The lattice community has opened the door for a very active collaboration between the theorists and experimentalists. With the EoS as an input, one can study the phenomenology of the hot and dense matter created at the heavy ion colliders. On the hand, there is a proposal of non-perturbative determination of the freezeout curve using lattice techniques, once the experimental data on cumulants of the charged hadrons are available.
A good understanding of the QCD phase diagram at zero baryon density has been achieved from the lattice studies. While the early universe transition from the QGP to the hadron phase is now known to be an analytic crossover and FIG. 18: The EoS for different isentropes using asqtad quarks are shown in the left panel, from [57]. In the right panel, the data for the energy density and pressure is compared for different isentropes using stout smeared staggered quarks, from [128]. not a real phase transition, it is observed that the chiral dynamics will have observable effects in the crossover region. One of the remnant effects of the chiral symmetry would be the presence of a critical end-point. The search for the still elusive critical endpoint is one of the focus areas of lattice studies, and the important developments made so far in this area are reviewed.
While QCD at small baryon density is reasonably well understood with lattice techniques, the physics of baryon rich systems cannot be formulated satisfactorily on the lattice due to the infamous sign-problem. A lot of conceptual work, in understanding the severity and consequences of the sign problem as well algorithmic developments in circumventing this problem is ongoing and is one of the challenging problems in the field of lattice thermodynamics.