Recent Progress on the QCD Phase Diagram

TL;DR

This paper reviews recent advances in understanding the QCD phase diagram through lattice QCD, highlighting progress in precision calculations of thermodynamic properties, topological fluctuations, and symmetry effects.

Contribution

It presents new insights from lattice QCD into the microscopic degrees of freedom, topological fluctuations, and the impact of anomalous symmetries on QCD phase transitions.

Findings

Enhanced precision in lattice QCD thermodynamics

Insights into topological fluctuations and screening masses

Understanding of U_A(1) symmetry effects on chiral transition

Abstract

Recent progress and the latest results on the bulk thermodynamic properties of QCD matter from lattice are reviewed. In particular, I will stress upon the fact that lattice techniques are now entering into precision era where they can provide us with new insights on even the microscopic degrees of freedom in different phases of QCD. I will discuss some instances, from the recent studies of topological fluctuations and screening masses. The progress towards understanding the effects of anomalous $U_A(1)$ symmetry on the chiral crossover transition and transport properties of QCD matter will also be discussed.

Full text

1 Introduction

The QCD phase diagram has driven the scientific curiosity of the community for more than thirty years, from its earliest versions discussed as early as 1980s. Understanding the phase diagram allows us to explain the origin of mass of 99.9% of the visible matter in the present universe. It has motivated large scale experiments from the LHC to the RHIC at BNL and is the special focus of the BES II runs at the RHIC during 2019-20. Several upcoming experiments at FAIR, NICA and JPARC are being designed to probe the phase diagram at very high baryon densities, yet to be understood. Experimental challenges aside, it is one of the most challenging problems in theoretical physics. Lattice studies have produced some of the remarkable results till now; it has now conclusively demonstrated that the phase transition at vanishingly small baryon densities is a smooth crossover [1, 2, 3]. Continuum results for bulk thermodynamic quantities like entropy density, pressure and the Equation of state (EoS) at zero baryon density are now known to very high precision [4] with new results on continuum estimates for the EoS available at baryon densities as large as $\mu_{B}/T\sim 2.5$ [5]. Efforts are underway to develop new lattice techniques to extend these calculations to even larger baryon densities $\mu_{B}/T\sim 3$. Moving ahead with these successes, I will show some instances, how lattice techniques are becoming mature enough to extend beyond bulk thermodynamic observables to understand the more microscopic details of the different phases of QCD, in particular, the microscopic origins of chiral symmetry breaking and deconfinement.

The review is organized as follows: In the first section, the recent updates on the thermodynamics crossover transition at $\mu_{B}=0$ are discussed. The anomalous $U(1)$ part of the softly broken chiral symmetry in QCD is believed to play an important role in determining the nature of the chiral phase transition in the limit when up and down quark masses are vanishingly small [6]. I will discuss the latest lattice results on the fate of $U_{A}(1)$ anomalous symmetry and how its origin of it can be traced back to the non-trivial topology of QCD. This leads to the next section which elaborates how the lattice community is trying to learn more about QCD phase diagram by varying the mass and the number of quark flavors within the so-called Columbia plot. The Columbia plot is now extensively studied along a new axis, by including an imaginary chemical potential to the QCD action. Finally I discuss how both imaginary chemical potential techniques as well as Taylor expansion in $\mu_{B}$ is allowing us to sketch the phase diagram in the finite density regime and possibly constrain a region in $T$-$\mu_{B}$ plane which may have the critical end-point. Other interesting topics discussed in the finite temperature sessions not included in this review are QCD at finite magnetic fields [7], strong coupling QCD [8] and QCD thermodynamics at large N [9].

2 Symmetries and phase diagram at $\mu_{B}=0$

Since the up and the down quark masses are much lighter than the intrinsic scale of QCD i.e $m_{l}=m_{u,d}<<\Lambda_{QCD}$, the $U_{L}(2)\times U_{R}(2)$ symmetry of the action is very mildly broken. $U_{L}(2)\times U_{R}(2)$ is isomorphic to $SU(2)_{V}\times SU(2)_{A}\times U_{B}(1)\times U_{A}(1)$ and 2+1 flavor QCD, has to a very good approximation, a $SU(2)_{V}\times SU(2)_{A}\times U_{B}(1)$ symmetry which is broken to $SU(2)_{V}\times U_{B}(1)$ leading to chiral symmetry breaking. The anomalous $U_{A}(1)$ part is always broken due to quantum effects. Though the chiral symmetry is exact in the limit $m_{u},m_{d}\rightarrow 0$, however remnants of it exist in chiral observables. For example, it was discussed in this conference [11] that the temperature at which an inflection point exist for the subtracted chiral condensate is consistent with the one at which the chiral susceptibility or its disconnected part peaks. An unweighted average of all these temperatures have been calculated in the continuum limit which allows for a more precise determination of the pseudo-critical temperature $T_{c}=156.5\pm 1.5$ MeV [10]. In contrast to the earlier reported value of $T_{c}=154(9)$ MeV [2] in the continuum limit, the systematic errors have reduced significantly by more than $80\%$. For recent updates on the results of chiral observables and measurement of $T_{c}$ with twisted mass fermions see [12].

The $U_{A}(1)$ part is an anomalous symmetry, there is no corresponding order parameter. From renormalization group studies of model quantum field theories with same symmetries as QCD, it has been observed that the order of phase transition for 2 flavor QCD depends on whether $U_{A}(1)$ breaking effects survive or gets effectively restored at $T_{c}$ [6]. Further studies with epsilon expansion [13] and conformal bootstrap [14] have revealed a possibility of a first order or even a second order phase transition of $U_{L}(2)\times U_{R}(2)/U_{V}(2)$ universality if the $U_{A}(1)$ is effectively restored near $T_{c}$ in contrast to an $O(4)$ second order transition if it remains broken. The magnitude of effective breaking of $U_{A}(1)$ can only be answered non-perturbatively and lattice techniques have immensely contributed towards a more systematic understanding of this issue. In order to quantify the effects of $U_{A}(1)$ at $T_{c}$, it was suggested quite sometime back to look at the degeneracy of the integrated two-point correlation functions of iso-triplet pseudo-scalar and scalar mesons [15]. The integrated correlation functions can be written in terms of the eigenvalues $\lambda$ and density $\rho(\lambda)$ of the QCD Dirac operator as $\chi_{\pi}-\chi_{\delta}=\int d\lambda\frac{4m_{l}^{2}\rho(\lambda)}{(\lambda^{2}+m_{l}^{2})^{2}}$, hence the properties of the eigenvalue spectrum as a function of temperature tells us about the fate of the $U_{A}(1)$. One way to trivially realize $U_{A}(1)$ restoration along with the chiral symmetry is to have $\rho(\lambda\rightarrow 0)=0$. On the other hand if the eigenvalue density has non-analyticities in its infra-red spectrum like $m_{l}^{\alpha}\delta(\lambda),\alpha\in[0,2)$ then $\chi_{\pi}-\chi_{\delta}$ is non-zero even in the chiral limit [16]. Recent theoretical studies suggest it is important to look at higher order correlation functions in all these mesonic quantum number channels [17]. In the chiral limit, calculations show that $U_{A}(1)$ breaking effects are invisible in upto 6-point correlation functions in the scalar-pseudo-scalar channel if the eigenvalue density goes as $\rho(\lambda)\sim\lambda^{3}$ [17]. The main issues on the study of $U_{A}(1)$ addressed and reported in this conference are,

• •

If one studies the eigenvalue spectrum of QCD at the physical point how does it quantitatively change as one goes towards the chiral limit. Are these spectra very different?

• •

Status of the finite volume and finite cut-off effects that crucially affects these studies.

New results on $\chi_{\pi}-\chi_{\delta}$ in 2-flavor QCD were presented in this conference [18], summarized in right panel of Fig. 1. It is observed that as one approaches the chiral limit, the finite volume effects could be milder (right panel of Fig. 1). For the physical quark masses the $U_{A}(1)$ breaking is still finite on lattices of size $48^{3}\times 12$ which seems to decrease to zero in the limit $m_{l}\rightarrow 0$. It would be interesting to study in detail how this reweighting of domain wall configurations work at large volumes and towards the chiral limit. The other approach reported was to calculate the eigenvalues of QCD Dirac operator for 2+1 flavors by fixing the strange quark mass $m_{s}$ to its physical value and reducing the light quark masses towards the chiral limit. New results on the eigenvalue spectrum of overlap Dirac operator on gauge ensembles generated using Highly Improved Staggered Quark (HISQ) discretization reported in Ref. [19], shows that the analytic part of the infrared spectrum is quite robust. $\rho(\lambda)\sim\lambda$ at around $1.1~{}T_{c}$ even when the light quark masses are reduced from $m_{l}=m_{s}/20$ [20] to $m_{s}/40$. A small non-analytic peak for $\lambda\rightarrow 0$ observed in the eigenvalue spectrum has been suspected due to effects of partial quenching [21]. In order of verify that the HISQ eigenvalue spectrum has been measured with the same valence and sea quark operators on fine lattices $64^{3}\times 16$ just above $T_{c}$. The small non-analytic peak seems to appear as one approaches the continuum limit [22]. It will be interesting to check this with other fermion discretizations like domain wall fermions also in the continuum limit though this study will be computationally much more intensive. As evident from the left panel of Fig. 1 both analytic and non-analytic parts of $\rho(\lambda)$ contribute to $U_{A}(1)$ breaking ($\chi_{\pi}-\chi_{\delta}$ has been renormalized appropriately to ameliorate the effects of partial quenching) which seems to survive even for $m_{l}=m_{s}/40$ at temperatures upto $1.1~{}T_{c}$. This suggests $U_{A}(1)$ breaking survives towards the chiral limit [19].

Another observable that measures the topological fluctuations of QCD vacuum is the topological susceptibility $\chi_{t}$. In LATTICE 2017, an extensive discussion of results from different groups suggest that for $T>3~{}T_{c}$ the temperature dependence of $\chi_{t}$ is consistent with the expectations from dilute instanton gas approximation (DIGA) [23, 24, 25] whereas non-trivial temperature dependence is seen for $T_{c}<T<3~{}T_{c}$ [26, 23]. New results with twisted mass fermions in $2+1+1$ QCD also confirms this overall picture [12]. Though the temperature dependence of $\chi_{t}$ agrees quite well with DIGA for $T>3~{}T_{c}$, its magnitude has to be scaled by a factor of $\sim 9$ to match with the leading order semi-classical result at $T\sim 450$ MeV [25]. This is due to the fact that the semi-classical result includes the color screening function at LO which has a slow convergence with the coupling [23]. It was argued that the semi-classical expansion of instanton action may not be as uncontrolled at $T\gtrsim 1.5$ GeV [27]. It would nevertheless be important to measure $\chi_{t}$ for $T>1$ GeV on the lattice to observe this convergence. However it is assuring that in the context of axion mass estimation, the temperature dependence of $\chi_{t}$ plays the decisive factor [23], changing the scale factor from $15$ to unity only changes the axion mass by $20\%$.

Interesting algorithmic developments have been reported since LATTICE 2017 to measure $\chi_{t}$ to very high temperatures [28, 29]. Since topological tunnelings become rarer as one goes to higher temperatures, one has to sample a large number of configurations to measure $\chi_{t}$ making the problem computationally challenging. It has been shown that sampling the ensembles with a reweighting factor with coarse-grained definition of topological charge, reduces the probability to get stuck at one topological sector. New results on continuum extrapolated $\chi_{t}$ for pure gauge theory (in left panel of Fig. 2) at $\sim 4~{}T_{d}$ was reported [29], which were calculated with very moderate computational efforts. Reweighting techniques have been applied to QCD with stout fermion discretization and $\chi_{t}$ has been calculated [30] after performing very careful finite volume and continuum extrapolation at $T\sim 450$ MeV, see central panel of Fig. 2. The results are consistent with the earlier results of $\chi_{t}$ calculated using a different reweighting technique [24] performed along the temperature axis starting from a low temperature ensemble [25]. Whereas this finite temperature reweighting is expected to work well for pure gauge theory where the temperature dependence is along the expectations of DIGA beyond $T_{c}$, it is more non-trivial to extend this technique in full QCD where the $T$-dependence is more intricate than naive DIGA below $2.5~{}T_{c}$. It is assuring that new techniques [30] confirm earlier reported results. Several other algorithms in the context of quantum mechanics [28] are discussed which have potential to be applied to QCD, some techniques discussed earlier like metadynamics [31] requires more extensive application. The endeavor towards measuring rare topological fluctuations at high temperature QCD has motivated development of new lattice techniques which can be applied to a more general problem when one approaches the continuum limit, when irrespective of the temperature, the ensembles get stuck in one topological sector.

New studies on understanding the topology in QCD near the chiral crossover region have been reported in this conference. In fact higher moments of the free energy $F(\theta)$ are more sensitive to the microscopic topological objects [32]. It has been reported earlier that the fourth moment of $F(\theta)$ has a value that is different from DIGA in the range $T_{c}~{}<T~{}<2~{}T_{c}$ [26]. This naturally leads to the question: what explains such an observation? At finite temperature the eigenvalues of the Polyakov loop at spatial infinity or the holonomy characterizes the properties of the instantons. For trivial holonomy the finite action solution at non-zero temperatures or calorons have been known for quite sometime [33]. Towards the end of 90’s, the calorons with non-trivial holonomy were discovered [34]. In fact it was shown that calorons in $SU(N)$ gauge theory consists of $N$ dyons, which carry a fraction $1/N$ of the net topological charge. Additionally dyons carry both color electric and magnetic charges and they combine in a way that the calorons are charge neutral objects. Calorons with trivial holonomy cannot explain the mechanism of confinement in gauge theories, mean-field studies of dyon gas hints to the fact that they may be a key towards understanding confinement [35]. It is therefore important not only to establish the existence of such objects non-perturbatively in QCD but understand their interactions.

An new study has been reported in this conference [36] furthering the earlier studies on dyons  [37, 38]. QCD ensembles are generated during a Monte-Carlo evolution with (anti)-periodic boundary conditions imposed along the temporal direction for (fermions) gauge fields hence an isolated dyon cannot exist on the lattice. However zero-modes of the valence Dirac operator with a general boundary condition $\psi(\tau+\beta)=\rm{e}^{i\phi}\psi(\tau)$, such that the twist angle $\phi$ lies between the eigenvalues of Polyakov loop, will detect the dyon characterized by the difference between these eigenvalues. This technique has been used to detect and characterize the zero modes of the overlap operator with different boundary conditions on Möbius domain wall fermion sea ensembles at temperatures between $T_{c}<T<1.1~{}T_{c}$ [36]. In fact density profiles of the zero mode wavefunctions show a very good agreement with the analytic profiles of dyons and their characteristic fall-off at large distances have been suggested as one of the signatures to identify dyons on the lattice. From a detailed study of near-zero modes, the interactions between dyons have been inferred qualitatively [36]. These insights will eventually lead us to an understanding of the mechanism of deconfinement and the yet un-explained temperature variation of $\chi_{t}$ just above $T_{c}$. At higher temperatures, $T>2~{}T_{c}$, the holonomy is trivial but there may be localized fluctuations of the Polyakov loop value, which is conjectured to provide the ’disordered’ landscape to localize the bulk eigenfunctions of the QCD Dirac operator [39]. A new study of the localization properties of the overlap Dirac operator on $2+1+1$ twisted mass sea ensembles (with pion mass of $\sim 370$ MeV) has been reported in this conference [40]. It further provides support for the conjecture that local negative fluctuations of the Polyakov loop provides the disorder required for localization of bulk eigenmodes and also reports that a dilute instanton gas cannot support such a localization [40].

Updates on the quark mass and volume dependence of $\chi_{t}$ in 2 flavor QCD with overlap fermions have been reported in this conference [41]. This study seems to suggest that the $\chi_{t}$ does not vanish linearly as $m_{q}$ but may either go as $m_{q}^{2}$ or rather abruptly vanishes for quark masses smaller than a critical mass $\lesssim 10$ MeV on a lattice of volume $~{}(2.4~{}fm)^{3}$. When the volume is increased to $\sim(3.6~{}fm)^{3}$, the value of $\chi_{t}$ at $m_{q}\sim 10$ MeV increases to a non-zero value, whereas for even smaller masses it seems to be consistent with zero with larger errors. The gluonic definition of $\chi_{t}$, however gives a non-zero value even for masses $m_{q}<5$ MeV. It would be interesting to check if this difference in values of $\chi_{t}$ measured using gluonic and fermion methods as a function of quark mass is resolved as one goes to the infinite volume and continuum limits.

Calculating bulk thermodynamic quantities of QCD on the lattice has interesting developments in past couple of years, both in terms of new techniques and results. There are updates on the EoS using non-perturbatively $\mathcal{O}(a)$ improved Wilson fermions where gradient flow is employed to fix the renormalization of the energy momentum tensor [42] on the lattice [43]. Since the continuum limit is not yet achieved, the idea is to look at the plateau of relevant quantities as a function of $a^{2}/t$ where $a$ is the lattice spacing and $t$ being the flow time [43]. The systematics of taking $t\rightarrow 0$ before continuum limit is studied in numerically inexpensive $SU(2)$ and $SU(3)$ gauge theories [44] and reported in this conference [45]. At present results for entropy density and subtracted chiral condensates on lattices of size $a\sim 0.09$ fm, are consistent with improved versions of staggered fermions whereas the interaction measure has still large errors. Update on the measurement of chiral susceptibility has also been reported [46]. New applications of gradient flow to fix the renormalization of energy momentum tensor correlators in full QCD with $m_{\pi}/m_{\rho}\sim 0.6$ has been discussed [47]. The correlators of $T_{12}$ show a plateau-like behavior near large Euclidean times $\tau\sim N_{\tau}/2$ for different flow times at $T=232$ MeV whereas for diagonal $T_{ii}$ correlators the plateau is quite noisy. Using model ansatz for spectral functions, the shear and bulk viscosities have been measured, the latter with larger errors. At present the results using the HTL ansatz cannot be differentiated from the Breit-Wigner ansatz for $T>200$ MeV and further spectral reconstructions are being studied [47]. A new result of the jet quenching parameter for $SU(3)$ gauge theory was reported [48]. Precise measurements of these real-time coefficients [49] will ultimately tell us how perturbative is the QCD medium beyond $T_{c}$. This is an evolving area, where lattice techniques need further development and has a promising potential. For recent updates on measurements of other real-time quantities like photon and di-lepton rates see Ref. [50] and a talk in this conference [51].

The EoS of QCD has now been measured with HISQ fermions for temperatures upto $2$ GeV by carefully performing continuum extrapolation [52], results of which were discussed in this conference (left panel of Fig. 2). The results are consistent with expectations from 3-loop HTL perturbation theory (without the static magnetic contribution). In fact measurements of the screening correlators at finite temperature for mesonic excitations in QCD [53] reveal that though in vector and axial-vector channels, the convergence to their perturbative estimates is quick, the scalar-pseudo-scalar excitations have a very slow convergence towards the perturbative value. Larger symmetries $SU(2N_{f})$ of fermion charge seem to be visible through the degeneracy of screening correlators of vector $V_{x}$ and tensor $T_{t}$ excitations as reported in [54] for $T>2~{}T_{c}$, when $U_{A}(1)$ is approximately restored. Near the perturbative regime at $T\sim~{}5~{}T_{c}$, these symmetries are again observed to be broken explicitly. All these studies hint to the fact that QCD medium is still non-perturbative beyond $T_{c}$ and the elementary excitations of the plasma have far more intricate structures than just free quark and gluon-like quasi-particles. The production of strange degrees of freedom is one of the other proposed signatures of a non-perturbative quark-gluon plasma. The FASTSUM collaboration have reported [55] on the parity restoration in different strange baryon channels near $T_{c}$. Though the $S=1$ baryon parity partners becomes degenerate like the non-strange baryons immediately near $T_{c}$, for higher strangeness sectors the parity restoration seem to occur much slowly, at $T>T_{c}$.

3 Towards Understanding the Columbia plot

A deeper understanding of the phase diagram of QCD is obtained when one looks at a more fundamental problem: what is the fate of ’chiral’ phase transition when the masses of quark flavors are varied. In left panel of Fig. 3 the current status of the famous Columbia plot is summarized. QCD with physical quark masses lie in the crossover region extended for a range of $m_{u},m_{s}$. The upper right corner of the plot is much better understood since for quark masses infinitely large, it corresponds to $SU(3)$ gauge theory which has a first order transition. This first order region is separated from the crossover region by a $Z(2)$ second order line. The lower left corner is comparatively much less understood. From model QFTs with same symmetries as $N_{f}=3$ QCD, it is expected that a first order region exist which should again be separated from the crossover region by a second order $Z(2)$ line. In fact $Z(2)$ scaling studies of chiral susceptibilities along the diagonal with $m_{s}=m_{u,d}$ on $N_{\tau}=6$ lattices with HISQ fermions constrain the $Z(2)$ line to exist for pion masses $m_{\pi}<50$ MeV [56]. With clover improved Wilson fermions the corresponding critical pion mass is $m_{\pi}<170$ MeV but for rather coarse lattices at present [57]. However for both staggered as well as Wilson fermions, the first order region tends to shrink when the lattice spacings are made finer. In a very insightful report [58] it has been motivated that as one approaches the continuum limit, the first order region for $N_{f}=3$ or even $N_{f}=4$ will shrink even further. For more updates on the status of critical $m_{\pi}$ for $N_{f}=4$ QCD with Wilson fermions, see [59]. The other question that naturally arises in this context, is whether this first order region end at a tricritical point for $m_{u,d}=0$ and a finite $m_{s}$ or continues all the way to the $m_{s}\rightarrow\infty$ axis. Which of these two scenarios survive in the continuum limit may ultimately be related to the fate of $U_{A}(1)$, which is not yet conclusively known. Already with coarser lattices, the first order region seems to be quite tiny in the lower left corner of the plot. If indeed this first order region survive as a tiny strip parallel to the $m_{s}$ axis and continue to $m_{s}\rightarrow\infty$ i.e. $N_{f}=2$ axis, then it is expected that the corresponding $m_{u,d}$ is much smaller than physical quark masses. The arrows on the plot in Fig. 3 indicate the directions of some of the current lattice studies in this regard. Summarizing them,

• •

The green arrow shows the JLQCD approach for $N_{f}=2$ QCD, explained in the previous section [18]. With the current lattice volume $(2.4~{}fm)^{3}$ and spacing $a^{-1}=2.6$ GeV, the results seem to suggest that $U_{A}(1)$ is restored for $m_{u,d}\lesssim 5$ MeV and could well be in the first order region. However the expectations from $N_{f}=3$ QCD seems to suggest that in the continuum, the first order region, if it survives and continues from the lower left corner all the way to the $N_{f}=2$ axis will very narrow characterized by $m_{u,d}<<5$ MeV. It will be important to reconcile both these results in the continuum limit. A related work discussed in this conference was to extract $T_{c}$ from a reweighted spectral density of QCD and thus obtain the order of transition in $m_{u,d}\rightarrow 0$ limit [60].

• •

The blue line on the Columbia plot shows the other approach by the HotQCD collaboration [61], where $m_{s}$ is fixed to its physical value and $m_{u,d}$ successively reduced to check whether one approaches the $Z(2)$ line to the left or goes over to a $O(4)$ second order line. New results on chiral susceptibility $\chi_{M}$ for $N_{\tau}=8,12$ lattices with HISQ fermions discussed in the conference, suggest that the peak of $\chi_{M}$ decreases with volume ruling out first order phase transition for $M_{\pi}>80$ MeV. Scaling studies of the chiral condensate normalized by $\chi_{M}$ seems to rule out $Z(2)$ scaling for $M_{\pi}>55$ MeV, see right panel of Fig. 4.

• •

New studies on the eigenvalue distribution of $2+1$ flavor QCD with HISQ fermions discussed earlier [19] also follows along this blue arrow. It has a different motivation, to find out if $U_{A}(1)$ remains broken as light quark mass is successively reduced from its physical value. If indeed $U_{A}(1)$ is broken, the $Z(2)$ line will not exist and one will directly hit the $O(4)$ line when moving towards $m_{u,d}\rightarrow 0$. The eigenvalue densities as observed for 3 pion masses $M_{\pi}\sim 160,140,110$ MeV seem to support this latter scenario.

A more general approach has been discussed in this conference: to vary the $N_{f}$ as a continuous parameter and study the fate of the chiral phase transition [62]. The idea is to start with $N_{f}=3$ QCD with finite quark masses in the first order region and zoom in to the tricritical scaling regime to extract $N_{f}^{tric}$ such that $m_{q}\sim(N_{f}-N_{f}^{tric})^{5/2}$. For $N_{\tau}=4$ lattices the $N_{f}^{tric}<2$, which seems to suggest a first order transition for $N_{f}=2$ (see right panel of Fig. 3). These results are being further verified in the continuum limit.

3.1 Adding a new axis to the Columbia plot

It was suggested in Ref. [63] that adding a third axis in form of an imaginary quark chemical potential $i\mu_{q}$ can further impose constraints on the 2D conventional Columbia plot. The QCD partition function in presence of $i\mu_{q}$ is free from ’sign-problem’ and has symmetries $Z(\frac{\mu_{q}}{T})=Z(-\frac{\mu_{q}}{T})$ and $Z(\frac{\mu_{q}}{T})=Z(\frac{\mu_{q}}{T}+\frac{2ni}{3}\pi)$ for $n\in\mathcal{Z}$. The center symmetry is thus again a good symmetry even in presence of finite quark masses. The phase of the Polyakov loop is an observable in this case, which will identify the different $Z(3)$ sectors as $i\mu_{q}$ is varied. For the Roberge-Weiss (RW) points characterized by $\mu_{q}=(2n+1)\frac{i\pi T}{3}$ [64], there is a transition between adjacent center sectors, which is first order for high temperatures and a smooth crossover for lower temperatures. From continuity, the first order lines should end in a second order RW end-point. The interesting question is how the deconfinement and chiral transitions at $\mu_{q}=0$ connect to the RW point. For heavy quark masses the first-order deconfinement lines, from reflection symmetry of the partition function, is expected to meet at the RW point, which will be a triple point. For intermediate values of $m_{q}$, the crossover curve at $\mu_{q}=0$ may meet at the RW end-point, expected to be in $Z(2)$ universality class. The chiral limit is however more interesting. Numerical simulations, initially on $N_{\tau}=4$ lattices with staggered fermions have shown a first order RW transition for both $N_{f}=2,3$ [65, 63], likely to survive in the chiral limit [67], confirmed later in studies with Wilson fermions [66]. This scenario is summarized in the modified Columbia plot at $\mu_{q}=\mu_{B}/3=i\pi T/3$, shown in left panel of Fig. 4. It is expected that the $Z(2)$ second order transition at intermediate masses is separated from the first order regions for both $N_{f}=2$ and $N_{f}=3$ by tricritical points. Now, what are its consequences for the $N_{f}=2$ chiral transition at $\mu_{q}=0$? If the $N_{f}=2$ chiral transition at $i\mu_{q}=0$ is,

a) second order, then the first order RW transition will end in a tricritical point for $\mu_{q}^{2}<0$.

b) first order, then the first order RW transition would end in a tricritical point at $\mu_{q}^{2}>0$.

The first lattice study along this line [67] was performed with staggered fermions on $N_{\tau}=4$ for different lattice volumes with $N_{s}=8,12,16$. The $Z(2)$ second order line was estimated for finite quark masses and for different values of $i\mu_{q}$ from Binder cumulants. Subsequently the $\mu_{q}^{tric}$ was estimated by looking at the tricritical scaling for $m_{u,d}$ in the chiral limit. The tricritical point was found at $\mu_{q}^{2}=0.85(5)T^{2}$ which seemed to suggest that $N_{f}=2$ chiral transition at $\mu_{q}=0$ is first order atleast on coarser lattices [67]. Subsequently improved versions of staggered fermions are been used to reduce lattice cut-off effects, which play a decisive role in this study. The most recent high statistics studies are being performed for $2+1$ flavor QCD by keeping the $m_{s}$ fixed to its physical value and reducing the $m_{u,d}$ at $\mu_{q}=i\pi T/3$ along the blue line shown on the lower RW plane of Fig. 4. To summarize these results:

• •

Studies with stout-smeared staggered fermions have been performed for several lattice spacings $N_{\tau}=4-10$ with current state of the art being $a=0.1$ fm [68]. The light quark masses have been varied such that the lowest pseudo-Goldstone pion mass achieved in the numerical studies is $50$ MeV. The largest volume is $N_{\sigma}=32$ such that $M_{\pi}L>1$. From scaling studies of the Polyakov loop susceptibility, a first order RW transition is not observed for $M_{\pi}\geq 50$ MeV, which in the continuum limit would imply that the first order region, if it continues to the $\mu_{q}=0$ plane would be a very narrow strip parallel to the $m_{s}$ axis. The other question addressed in this study is how close are the RW and the chiral transitions. As evident from right panel of Fig. 5, the chiral and RW transition seem to follow each other as one reduces the $m_{u,d}$. The scaling studies of the subtracted chiral condensate near the RW point at present cannot distinguish between $O(2)$ universality scenario for $N_{f}=2$ and the $Z(2)$ universality expected at the RW transition [68].

• •

The RW transition is related to the restoration of $Z(2)$ symmetry. Under $Z(2)$ transformation, the real part of Polyakov loop does not change sign whereas its imaginary part changes sign. Hence the expectation value $\langle|ImL|\rangle$ is a good order parameter and will show $Z(2)$ scaling. Scaling studies performed with HISQ fermion discretization for $N_{\tau}=4$ and $N_{\sigma}=8-24$ around the chiral crossover transition temperature $T_{c}\sim 200$ MeV for $M_{\pi}=135-90$ MeV has been reported in this conference [69]. As evident from left and central panels of Fig. 5 a beautiful agreement with second order $Z(2)$ scaling is observed both for the order parameter and its susceptibility again confirming the previous independent finding [68].

The results using different improved versions of staggered fermions are converging to an agreement with no indication for a first order transition in the vicinity of RW fixed points for $M_{\pi}\gtrsim 50$ MeV. However as argued in [68], the other pion states in these studies are still quite heavy so it is important to revisit these studies on more finer lattices or with other fermion discretizations.

4 Status of QCD Phase diagram at finite density

Simulating QCD at finite density on the lattice is one of the most challenging problems in theoretical physics. New computational techniques and algorithms have been discussed in this conference in order to ultimately simulate dense and (cold) quark matter and understand the yet unexplored regions of the phase diagram. For relevant references and phenomenological applications for QCD at finite density, see the plenary talk by C. Ratti in this conference [70].

If indeed a first order transition occur in cold and dense QCD following clues from Nambu-Jona-Lasinio model, it should end in a critical end-point since we now know for sure that there is a crossover transition at $\mu_{B}=0$. Lattice is essential to establish if a critical end-point exist and to draw lines separating the phases of dense QCD matter. In this section, I will rather discuss how existing lattice techniques are allowing us to draw the chiral crossover line at small $\mu_{B}$ and what promise it holds to reach all the way to the critical end-point. Out of the many methods developed over the years to circumvent the sign problem, two of them have now been adapted for simulations at large volumes and towards the continuum. One of them is to simulate QCD at imaginary $\mu_{B}<\mu_{B}^{RW}$, calculate thermodynamic quantities like baryon number density and extrapolate to the real $\mu_{B}$ plane [71, 72]. The other method is to calculate the partition function at $\mu_{B}=0$ and expand it as a Taylor series in $\mu_{B}$ [73]. If indeed singularities like a critical end-point exist in the $T-\mu_{B}$ plane, then its location will determine the radius of convergence of the series [74].

Both methods have been used to calculate the curvature $\kappa_{2}$ and higher derivatives $\kappa_{4}$ of the chiral crossover line at small $\mu_{B}$, defined through $\frac{T_{c}(\mu_{B})}{T_{c}(0)}=1-\kappa_{2}\frac{\mu_{B}^{2}}{T_{c}(0)^{2}}-\kappa_{4}\frac{\mu_{B}^{4}}{T_{c}(0)^{4}}$. The results using Taylor expansion of chiral condensate with HISQ fermions for different $\mu_{X}$ where $X$ represents quantum numbers like baryon no., strangeness etc. were discussed in this conference [11] and summarized in the left panel of Fig. 6. The status of all recent lattice studies is summarized succinctly in the right panel of Fig. 6 from Quark Matter 18 review by M. D’Elia [75]. For quite a few years, there was an apparent disagreement between the values of $\kappa_{2}^{B}$ obtained using Taylor expansion and imaginary $\mu$ methods. A careful continuum extrapolation was the key to this resolution [76]; the new continuum results with HISQ fermions uses $N_{\tau}=6,8,12,16$ [10] data and with stout-smeared staggered fermions [76] uses $N_{\tau}=6,8,10$ results. The values of $\kappa_{2}^{B}\sim 0.01$ and $\kappa_{4}^{B}\sim 0$ suggest that the pseudo-critical line is almost flat for small $\mu_{B}$ and bends inwards very gradually towards larger $\mu_{B}$. Observables like the chiral disconnected susceptibility shown in left panel of Fig. 7 from Ref. [11, 10] also show a very mild dependence on $\mu_{B}$ for $\mu_{B}<250$ MeV.

New results for higher order fluctuations of conserved quantum numbers: baryon, charge, strangeness (B,Q,S) calculated using both these methods are available. For the imaginary $\mu$ method the latest high-statistics results are available from two different groups, both in $2+1$ QCD with stout-smeared staggered quarks. In one of these studies [77], all possible diagonal and off-diagonal second order susceptibilities in the $(i\mu_{B},i\mu_{Q},i\mu_{S})$ plane were calculated on $32^{3}\times 8$ lattice for temperatures between 135-350 MeV. Approximating these second order correlations and fluctuations by a polynomial of $\mathcal{O}(\mu_{B}^{i}\mu_{S}^{j}\mu_{Q}^{k})~{},~{}i+j+k\leq 8$ and extrapolating to the real plane, all higher order susceptibilities upto 8th order have been calculated. The other group uses a finer $48^{3}\times 12$ lattice to calculate all possible correlation and fluctuations of B, Q, S for temperatures between 135-220 MeV at 8 different imaginary $\mu$ values [78]. These were fitted to a polynomial of $\mathcal{O}(\mu_{B}^{i}\mu_{S}^{j}\mu_{Q}^{k})~{},~{}i+j+k\leq 10$, where the eighth and tenth order data were put in as priors. Using this fitting procedure, higher order susceptibilities upto 8th order were reported. Since there is a discontinuity of the imaginary baryon number density at the first RW point for $T\geq T_{RW}$, this naturally limits the number of imaginary $\mu$’s where the simulations can be performed and hence the extrapolation to real $\mu$. Given that range of simulations are limited to $\mu_{B}/T\in[0,i\pi)$, this method works better for $T<T_{c}$ but the systematic errors start dominating for $T>T_{c}$. In the Taylor expansion method, the pressure is expanded as a series in $\mu_{B}/T$ where the expansion coefficients are $\mu_{B}$-derivatives of pressure i.e. the higher order susceptibilities, calculated at $\mu_{B}=0$. These quantities involve derivatives of Dirac operator and each derivative is associated with an inverse of the Dirac matrix. The higher order fluctuations thus contains many such terms with alternating signs for subtle cancellations of the divergences to give a finite result. For $\chi_{6}^{B}$ or higher, the divergences may not exist which allows using a different technique [79] to compute them, which is computationally much cheaper compared to the conventional method [80]. A new numerical implementation of this technique which may allow to calculate even higher order fluctuations was discussed [81]. The current state-of-the-art results using Taylor expansion are correlations and fluctuations upto sixth order using Highly Improved Staggered quarks [5] and upto eighth order using unimproved staggered fermions [82]. Even with susceptibilities upto $\mathcal{O}(\mu_{B}^{6})$, the QCD EoS show very good convergence for $\mu_{B}/T\lesssim 2.5$, the continuum estimates of which can be found in [5].

For locating the critical end-point, the radius of convergence (RC) of the Taylor series of pressure or the baryon number fluctuation, $\chi_{2}^{B}$ has to be estimated. From the definition of RC $r_{2n}\equiv\sqrt{2n(2n-1)\left|\frac{\chi_{2n}^{B}}{\chi_{2n+2}^{B}}\right|}$ it is not known a-priori how large the order $n$ should be chosen in order to reliably extract this quantity on the lattice. The current estimates of the radius of convergence is summarized in right panel of Fig. 7. Most of them except the reweighting data are from $N_{\tau}=8$ lattices. The $r_{4}$ already deviates from the Hadron Resonance gas model (HRG) estimates by $\sim 30\%$ for $T\sim 145$ MeV [5]. There is a substantial difference between the $r_{2}$ and $r_{4}$ estimates so one needs atleast $r_{6}$ to get a reliable prediction for the RC. The yellow ’exclusion’ region comes from the upper error bar of the $\chi_{6}^{B}$ measured using the HISQ fermions whose central values are given by the blue points [5, 83]. Results using stout-smeared staggered quarks [77] also favor a larger $\mu_{B}^{CEP}/T$ than using the standard staggered quarks [82] shown by the black solid point. All these results should ultimately agree in the continuum limit. The $\mu_{B}^{CEP}/T$ using reweighting techniques from Ref. [84] favors a lower value; it will be interesting to confirm this in the thermodynamic limit. To summarize, the present lattice data for $\chi_{n}^{B}$ already deviates from naive expectations from HRG model at $T>145$ MeV; the higher the order $n$, the more visible is the deviation. Moreover the present lattice data favor a small curvature of the chiral crossover line [76, 11]. Furthermore the fact that $\kappa_{4}^{B}\sim 0$ suggest if a CEP exist in the phase diagram then $T_{CEP}\lesssim T_{c}(\mu_{B}=0)$. For the case $T_{CEP}/T_{c}(0)\sim 0.95$, lattice data already suggests stronger departure from HRG results and can at present provide a suggestive lower bound, $\mu_{B}^{CEP}>4T$ [83]; the convergence to the actual value will depend on a more precise calculation of $\chi^{B}_{8}$. On the other hand if $\kappa_{6}^{B}$ or $\kappa_{8}^{B}$ have strong contribution to the curvature of the pseudo-critical line and $T_{CEP}/T_{c}(0)\leq 0.9$, the RC estimates would be more closer to HRG values, hence will show extremely slow convergence as a function of the order $n$.

5 Outlook

In this review, I hopefully could convince that the quest to understand the phase diagram of QCD has led to the many interesting theoretical and algorithmic developments in lattice gauge theory. Efforts to understand the chiral transition better, has led to a rich theoretical knowledge of QCD in the $m_{u,d}$-$m_{s}$ plane, additionally along the imaginary chemical potential as a new axis and now even as a function of $N_{f}$ as a continuous parameter. The role of $U_{A}(1)$ anomaly on the chiral phase transition is not yet fully understood but has led to many new insights on the microscopics of the QCD Dirac operator and its intimate connections to topology. The topological structures in QCD and their interactions is long suspected to drive chiral symmetry breaking and confinement at finite temperature and/or densities; new insights on which are coming from lattice studies. Moreover a strong motivation to quantify topological fluctuations at high temperatures have led to recent developments of interesting new algorithms that has even more wider applicability for lattice simulations near the continuum limit. The quest to go deeper along the $\mu_{B}$ axis of the phase diagram has triggered developments of algorithms and new techniques to circumvent the sign problem, with initial bounds available from lattice towards constraining the location of the critical end-point. The EoS characterizing different phases of QCD upto $\mu_{B}/T\lesssim 2.5$ is now available; an increasing sophistication of lattice techniques is leading towards quantifying its dynamical properties.

Acknowledgements I thank the organizers for the kind invitation. Also I express my gratitude to Swagato Mukherjee who kindly agreed to fill in for my absence at the conference due to unavoidable circumstances, at a short notice and share his perspectives. I thank all the researchers working in finite temperature and density lattice QCD who have generously shared their inputs and results, specially to Massimo D’Elia, and also for the exciting talks in LATTICE 18; it indeed shows how vibrant this field of research is. Last but not the least, I am grateful to the members of HotQCD collaboration for a very enjoyable collaboration over these years and their helpful feedback.