Frontiers of finite temperature lattice QCD

I review a selection of recent finite temperature lattice results of the past years. First I discuss the extension of the equation of state towards high temperatures and fi- nite densities, then I show recent results on the QCD topological susceptibility at high temperatures and highlight its relevance for dark matter search.


Introduction
Quantum chromodynamics (QCD) is the underlying quantum field theory of the strong interaction that is probed in heavy ion experiments, such as the Relativistic Heavy Ion Collider (RHIC) at Brookhaven, and at the Large Hadron Collider (LHC) at CERN. The relevance of QCD predictions, however, stretches out to temperatures much beyond the reach of these experiments, and must be considered in cosmological scenarios of the Early Universe [1].
Lattice simulations provide an excellent tool to study equilibrium QCD at all but the perturbative temperature ranges of interest. The algorithmic difficulties that prevented the use of the light physical quark masses have been mostly overcome [2], and the systematics of the lattice discretization is under control by using sufficiently fine lattices and performing a continuum extrapolation, allowing lattice QCD to reconstruct the hadron spectrum [3,4].
Thermodynamics calculations are also possible with physical quark masses with either staggered [5] or domain wall fermions [6]. The transition is a cross-over [7], its temperature has been calculated based on the chiral observables, T c ≈ 155 MeV [8][9][10][11], the deconfinement aspects did not allow a clear-cut T c definition, though all possible definitions are within the peak range of the chiral susceptibility [10]. Recently the single-quark entropy was introduced based on the Polyakov loop which also peaks near the chiral T c [12]. The equation of state has been determined well beyond the transition regime, including the temperature derivatives, such as the trace anomaly and higher derivatives, like the speed of sound and heat capacity [13][14][15]. In quenched QCD the temperature range stretches out into the perturbative regime [16,17].
Several experimental as well as lattice papers have been lately devoted to the fluctuations of conserved charges [18,19]. In a grand canonical ensemble (with zero or non-zero chemical potentials) conserved quantum numbers, such as the baryon number, electric charge or strangeness, fluctuate according to equilibrium thermodynamic features, that are given by even derivatives of the free energy with respect to the corresponding chemical potentials. The odd derivatives vanish due to the C-symmetry of the QCD action. The magnitude of these fluctuations are distinctive features of the a e-mail: borsanyi@uni-wupperta.l.de arXiv:1612.06755v1 [hep-lat] 20 Dec 2016 EPJ Web of Conferences two phases around the QCD transition [20]. Precision experimental studies provide the values of ratios of such fluctuations for several collision energies as part of the RHIC beam energy scan program [21,22]. The obtained data can be described by a grand canonical ensemble assuming a temperature and chemical potential which are referred to as chemical freeze-out parameters. This is a generalization of the earlier statistical hadronization model, which was based on particle yields instead of fluctuations [23,24]. The fluctuation-based freeze-out temperatures lie in general below the original yield-based freeze-out curve [25].
These achievements of lattice QCD are mostly limited to zero chemical potential and to neartransition temperatures. In this work we describe the progress towards finite densities, relevant for RHIC phenomenology, and towards higher temperatures which links QCD into the cosmological context. Of particular recent interest is the lattice calculation of the topological susceptibility in QCD, which provides the non-perturbative theoretical input for axion cosmology.

Equation of state at high temperatures
Even in the absence of quarks the study of high temperature lattice fields becomes technically very difficult if the simulation volume is not to shrink with the temperature. Typical simulation programs use a given a lattice geometry of N τ temporal slices in the Euclidean time direction and N σ sites in each spatial directions. On an isotropic lattice with spacing a, the temperature is T = 1/(N τ a), and the box length is given by L = aN σ . Thus, with fixed geometry L drops with 1/T as L = N σ /(N τ T ). In the pure SU(3) theory Refs. [16,17] used very large spatial lattices so that the simulations could be carried beyond 10T c . While Ref. [16] used the conventional technique [26] where the trace anomaly is calculated from a pair of simulations at each temperature, the approach of Ref. [17] defines the ensemble in moving frame, allowing the extraction of the entropy from the off-diagonal elements of the energy-momentum tensor [27]. The latter technique requires an elaborate renormalization procedure, which was worked out in Ref. [28]. From the entropy function the other observables are found through thermodynamic identities. We show the results in Fig. 1. Refs. [16,17] agree for high temperatures. Below 3T c disagreements are reported on the 4% percent level [29]. Lacking systematic error budgets the significance of these are difficult to assess.
A further interesting method was used in Ref. [30], which is based on the gradient flow. It was observed the signal-to-noise ratio in the lattice energy-momentum operator can be significantly enhanced by a filtering of the UV-modes, resulting in automatically renormalized observables in the gauge sector [31,32]. A finite energy-momentum tensor was obtained at positive flow times perturbatively [33], these findings have already been generalized to full QCD [34]. The first numerical results with Wilson quarks have been published in Ref. [35]. Though these directions are very promising, an implementation with physical quarks masses is still missing. So that we can discuss continuum extrapolated results with a physical particle spectrum we turn to the staggered formulation.
The study of the equation of state in full QCD is by far more challenging. In Ref. [13] a continuum estimate was given using stout-staggered quarks with physical quark masses. Later almost the same results were obtained after continuum extrapolation [14] or with a different type of staggered quarks [15]. These results are summarized in Fig. 2.
The studies presented in Fig. 2 were based on the simulations of two discretizations of the 2+1 flavor action. Thus the effect of the charm was here neglected, although from about 300 MeV temperature a significant contribution is expected based on perturbative results [36]. Earlier partially quenched works suggested an earlier onset of the charm [13,37], which was later explained as a renormalization effect [38]: in the equation of state determination a key role is played by the a(β) scale function (with the inverse coupling β = 6/g 2 0 ), which, too, is significantly modified by the inclusion of the charm quark.  [16], this plot uses the non-standard normalization ( − 3p)/T 2 . The entropy density (blue), pressure (green) and the energy density (red) on the right panel is taken from Ref. [17], which is the latest publication on the subject.  Continuum results for the pressure (red), energy density (green) and trace anomaly (blue) as a function of temperature, as published by the Wuppertal-Budapest group [14] using the 2stout action and the HotQCD collaboration [15] using the HISQ action. Right: The speed of sound (c 2 s = d p/dε) plotted as a function of the energy density. On the top axis we label the respective temperatures.
A satisfactory treatment of the charm degree of freedom, especially if fine lattices are desired, required the introduction of a new lattice action with dynamical charm. In Ref. [38] the already used 2-stout staggered action was equipped with a charm quark. A year later a new large scale staggered program was launched by the Wuppertal-Budapest group, using four steps of stout smearing [39], the preliminary results for the equation of state were shown at the Quark Matter conference in 2012 [40]. The tuning procedure was refined and the finite temperature ensembles were introduced in Ref. [41]. The ETM collaboration also introduced the charm quark as the counterpart of the strange quark, covering a pion mass range between 210 and 470 MeV. They calculated the equation of state with 2+1+1 dynamical flavours for 370 MeV pions and three lattice spacings [42]. The MILC collaboration used the HISQ action with 2+1+1 dynamical flavors and calculated the trace anomaly in Ref. [43]. They reported very small discretization effects with lattices as coarse as N τ = 6, 8 and 10, the data were taken with a pion mass of 320 MeV. This effort was, however, not pursued further towards the physical point.
In the rest of this section we present the result of the Wuppertal-Budapest group, that published continuum extrapolated results with physical quark masses in Ref. [44]. For the equation of state this paper used the same scheme as earlier in the quenched case [16]. Each finite temperature ensemble is accompanied with an other ensemble with the same bare parameters, but belonging to a different temperature (different N τ ). The quartic divergent contributions are equal in the two simulations, the difference in the trace anomaly (I = −3p) contributions is [I(T )−I(T/2)]/T 4 assuming a factor of two in N τ between the two ensembles. These ensemble pairs come supplementary to the [I(T ) − I(0)]/T 4 pairs, which are conventionally used, since at zero temperature the trace anomaly vanishes I(0) = 0. The spatial box size L has been kept large even with increasing temperatures, such that LT c 2. E.g. at T = 880 MeV we used the geometries 64 3 × 6, 96 3 × 8, 128 3 × 10 and 128 3 × 12, and with the same bare parameters I(T/2) was calculated on 64 3 × 12, 96 3 × 16, 128 3 × 20 and 128 3 × 24 lattices. The ensembles up to 600 MeV temperature were listed in Ref. [41]. The to find the absolute I(T ) result where n is the smallest non-negative integer with T/2 n < 250 MeV. A comparison to the 2013 data set of Ref. [14] shows that the effect of the charm is increasingly relevant from about 300 MeV temperature.   It in natural to expect that towards 1 GeV temperature perturbation theory becomes applicable. The low orders of the conventional perturbative series fail to show convergence. Up to O(g 6 log g) the QCD pressure is known analytically [45], the O(g 6 ) term can be found numerically through dimensional reduction to the three dimensional effective theory [46]. Here the idea of fitting the O(g 6 ) term directly to 4D data was also discussed for the quark-less theory. This direct fit was later recalculated as high temperature data became available [16]. A similar strategy was implemented now in the CONF12 unquenched case. However, because of the non-negligible charm quark mass the result could not be fitted by the four-flavor perturbative result. Following a tree-level version of the ratio in Ref. [36] one writes Here SB (3) and SB(4) are the Stefan-Boltzmann limits of the three and four flavour theories, respectively, and F Q (m c )/T gives the free energy for one free charm quark with m c = 1.29 GeV. The leading perturbative corrections to Eq. (2) are calculated in Ref. [36], these are well below the percent level near the highest simulation point at about 1 GeV. Thus we attempt to give a perturbative description to our simulation results by fitting Eq. (2) to the pressure result at a single temperature T max = 1 GeV. Here we use for p (2+1+1) (T ) m c =0 the known O(g 6 log g) formula plus one fit parameter (q c ) controlling the g 6 contribution. The coupling is translated to temperature using the four-loop running of the coupling with Λ MS = 290 MeV, the standard value for n f = 4 [47].  The non-trivial result is shown in Fig. 4. Fixing the only fit parameter to q c = −3000 (a parameter not far from that in the pure Yang-Mills theory in Ref. [16]), we have good agreement for both the pressure and the trace anomaly starting already from approx 500 MeV. There is actually a ±10% window in q c where a simultaneous agreement with the trace anomaly and pressure is granted. Since at high temperature the fitted formula is asymptotically equivalent to the known perturbative result, one may expect that the agreement persists for higher temperatures, too. A direct cross-check to lattice data beyond 1 GeV is, however, technically demanding. Although deep in the deconfined phase the simulation with dynamical quarks are inexpensive, the autocorrelation times are increasingly difficult to control. In addition, the systematics of the tuning of the bare parameters (or in other words, the setting of the scale and the quark masses) cannot be feasibly backed by zero-temperature simulations.
But once we demonstrated that simple parametrizations like Eq. (2) can be used to describe the quark mass effects in the deeply deconfined phase we can attempt to go for higher temperatures and add the bottom quark using the following ratio where m b (m b ) = 4.18 GeV is the mass of the bottom quark [47]. Eq. (3) is expected to work beyond the bottom threshold. Actually, the ratio of the four and five flavour perturbative results are equal the respective ratio of their Stefan-Boltzmann limit to about 0.3% accuracy in the relevant temperature range. The blue line in Fig. 4 represent the l.h.s. of Eq. (3). This is our best estimate to the QCD equation of state from 500 MeV temperature. More details of the fitting procedure are given in Ref. [44].
In the same paper the pressure and energy density contributions of photons, neutrinos, charged leptons are also added up. The electroweak sector, including the W ± , Z 0 and Higgs bosons as well as the top quark are taken from Ref. [48]. The full equation of state can then be used to close the Friedmann equations and to determine the cooling rate of the Universe:

Equation of state at finite density
An other extension of previous equation of state simulations is towards finite densities. RHIC events freeze-out in a chemical potential range between 20 and 400 MeV [24,49]. Thus for the phenomenology of the heavy ion collisions the exploration of the range up to µ B /T 3 is sufficient, for quark chemical potential this means µ ud /T 1. Conventional Monte Carlo simulations at nonvanishing real chemical potentials are not possible, since then the formerly positive definite fermion determinant may take complex values, and thus the exponentialized action cannot be interpreted as a probability weight. Moving the complex part into the observable introduces heavily oscillating contributions, which is known as the sign problem. On top of that, the most relevant configurations are probably never sampled due to the ill-tuned action, this is the overlap problem. An appropriate redefinition of the degrees of freedom have solved the complex action problem in many other theories [50], and recently much progress has been made to solve the overlap problem in gauge theories using density of states methods [51]. Other direct methods to solve the complex action methods include Lefschetz thimbles [52] and complex Langevin simulations [53,54]. Promising these new directions may be, their applicability to full QCD has not yet been proven.
For small chemical potentials, however, the leading coefficients of a Taylor expansion may already provide a satisfactory description of this range of the QCD phase diagram. Calculations to µ 6 order have already been available for a decade [55], the leading µ 2 order was explored two decades ago [56]. These works were neither with physical parameters, nor continuum extrapolated, though. It turned out that the technical difficulty of getting the Taylor coefficients of the free energy greatly increased as one moved on to lighter pions and to finer lattices. The reduction of the lattice spacing and the introduction of fat links (HISQ or stout) have finally allowed the calculation of the continuum limits of the leading µ derivatives of the effective action [57,58] as well as the leading Taylor coefficient for the equation of state [59].
It has been a great challenge to achieve an acceptable signal for the higher coefficients in realistic lattice simulations. For this the fourth and sixth order fluctuations of conserved charges had to be calculated. The determination of these fluctuations have been also motivated by the experimental availability of ratios of cumulants [21,22]. The first applications of lattice QCD have provided estimates for the freeze-out parameters [60][61][62]. The extrapolation of the fluctuations to finite chemical potential may allow the extraction of the curvature of the freeze-out line from fluctuation data [63][64][65].
The fourth-order derivatives in the continuum limit have been published in Ref. [41]. In Fig. 5 we show χ U 4 , the fourth derivative of the free energy density with respect to the up quark chemical potential (the kurtosis of the up quark number) from Ref. [66] and the analogous quantity for the CONF12 baryon number χ B 4 from Ref. [41]. In both cases the data are compared to the predictions based on hard thermal loop perturbation theory [67] and to those from dimensional reduction [68].  [66] and the analogous fluctuation for the baryon number (right) [41]. In both plots the data are confronted to the expectations based on HTL perturbation theory [67] and dimensional reduction [68].

T [MeV]
From the error bars in Fig. 5 one can quickly notice that the calculation costs of the fluctuations (and, hence, for the Taylor coefficients, too) rapidly increases as the temperature is lowered towards and below T c . To meet the challenge, the BNL-Bielefeld-CCNU group invested in a large pool (∼ 10 5 ) of configurations for each temperature, although mostly at the resolution N τ = 8. The Wuppertal-Budapest group collected statistics at N τ = 10, 12 and 16, mainly using imaginary chemical potentials in the simulation parameters. There is no complex action problem hindering simulations with imaginary values of the chemical potentials. As it was recently emphasized in Ref. [69] knowing the µ B -dependence of low order µ B derivatives of the free energy allows the calculation of higher derivatives. This principle was exploited by the Wuppertal-Budapest group in Ref. [70]. Both in the Wuppertal-Budapest simulation program [64,71] as well as in the BNL-Bielefeld-CCNU effort [65,72] care has been taken to correctly take into account the mixing of the strange and baryon chemical potentials. To reproduce the experimental setting, the expectation value of the strangeness has to be zero, and the electric charge is set to 0.4 times the baryon number to reflect the flavor imbalance of the projectile ions. In Ref. [70] it meant to match the imaginary strangeness chemical potential to the imaginary baryon chemical potential, thus setting a trajectory in the imaginary µ B − µ S plane. The Taylor coefficients then correspond to the total derivatives along this trajectory. The results are plotted in Fig. 6.

Topological susceptibility at high temperatures
Perhaps the most exciting topic of the past year was the topological susceptibility of QCD at finite temperature. It is defined as the second derivative of the free energy with respect to the θ parameter, which is introduced in the Euclidean theory as A nonvanishing θ parameter would introduce an electric dipole moment for the neutron [73,74], but this is constrained experimentally [75], such that |θ| < 10 −9 . Given the CP-violation coming from the CKM-matrix a fine-tuning is necessary to maintain a near-zero θ term, unless at least one quark has  Figure 6. Taylor coefficients of the QCD pressure in the transition region calculated from imaginary chemical potential data by the Wuppertal-Budapest group [70]. Similar results were presented at the Confinement12 conference by C. Schmidt using BNL-Bielefeld-CCNU data at µ B = 0. zero mass, in which case the θ term can be removed by chiral rotation. However, there is no quark with zero mass [76].
To solve this fine-tuning problem, the strong CP problem, Peccei and Quinn have proposed to introduce dynamics to θ [77]. Assuming a spontaneously broken global U(1) symmetry for a new dynamical field, the angular degree of freedom will the play of θ. This pseudo-Goldstone boson acquires a small mass (m a ) by the QCD chiral anomaly. This mass is then proportional to the second derivative of the effective potential: where f a is the U(1) symmetry breaking scale, yet to be determined. This angular degree of freedom is the axion field (for a brief review see Ref. [78]). The topological susceptibility χ t at zero temperature is well predicted by chiral perturbation theory [79], and can be calculated in quenched [80,81] and in the unquenched theory [82,83]. Early lattice studies have shown that the topological susceptibility sharply drops at the transition temperature [84,85]. The strong temperature dependence can be explained by the dilute instanton gas approximation (DIGA) [86,87]. In a simplified picture the action of a single instanton 2π/α s contributes as to the topological susceptibility. From the one-loop running of the coupling we have Actually, in the presence of fermions the configurations with non-trivial topology are suppressed by a factor of the light quark mass, which, at high temperature, is made dimensionless as ∼ m/T . The contributions together give a power law: For a full calculation see Ref. [88].
To demonstrate the power law behaviour on lattice was a challenge even in the quenched theory [89][90][91]. The result is shown on the left panel of Fig. 8. The power law χ ∼ T −7 behaviour means that as the temperature is increased in a constant volume the weight of the configurations with nonzero topology shrinks with the same power, and the length of the Monte Carlo simulations must be extended accordingly. Moreover, an algorithmic difficulty, the topological freezing, preventing the crossing between topological sectors, worsens the problem. Suggestions include artificial decreasing of the weight of the trivial sector [90] and the use of the topological charge density correlator [92]. In Refs. [44,93] the following new strategy was introduced. The partition sum is actually a sum over the contributions from different topological sectors with charge Q: In a given volume at high temperature Z 0 Z 1 Z 2 . . . , and through P-symmetry: Z 1 = Z −1 , Z 2 = Z −2 , . . . If we keep the simulation aspect ratio R = LT fixed, the susceptibility can be written as Thus the temperature dependence of the topological susceptibility can be deduced from the trace anomaly difference between the sectors.
This way the exponent of the power law behaviour is obtained directly. In full QCD this amounts to: Here we used the notation · 1−0 for the difference between the expectation values in the Q = 1 and Q = 0 sectors. However, the generalization to full QCD is not as straightforward. Discretization errors are severe, they can only be slightly reduced by advantageous normalization, e.g. to the zero temperature EPJ Web of Conferences  Here a pre-factor had to be introduced for DIGA to match the lattice data [91]. Right: A compilation of recent results for full QCD. Magenta band: Bonati et al [94] using 2-stout staggered fermions and the gluonic definition of the topological charge. Blue and green bands by Petreczky et al [95] come from the gluonic definition and the disconnected chiral condensate, respectively, using HISQ fermions. The red band by the Wuppertal-Budapest group [44] comes from a combination of simulations with staggered and overlap quarks. χ t (T = 0) at the same lattice spacing [44,94]. One very important lattice artefact in the staggered formalism is the enhancement of the near-zero eigenvalues in the Dirac operator. The number of the zero modes is related to the topological charge though the index theorem. Even a large discretization error on a few eigenvalues in the scarce non-trivial configurations has probably small effect on other bulk quantities, for the topological susceptibility, however, the very weight of these configurations matter. To correct the weight, a reweighting strategy was suggested in Ref. [44]. Depending on Q, the staggered eigenvalues of the appropriate number of low modes are replaced by the would-be eigenvalues, given by the mass. Since only the weight of the configurations are affected, it is irrelevant whether the would-be zero mode or a mode with a very similar eigenvalue was picked. Note that the fermionic observables, e.g. ψ ψ f 1−0 also need to be redefined according the modified action [44]. In Fig. 7 we show examples for the improved continuum scaling as well as a brief finite volume study.
The lattice results on the topological susceptibility are summarized in Fig. 8. The left panel shows the continuum extrapolated quenched result together with the DIGA prediction. Although DIGA did predict the exponent correctly, an order of magnitude mismatch is found in pre-factor [91]. The full QCD result by three collaborations [44,94,95] is shown in the right panel of Fig. 8.
Having determined the QCD topological susceptibility over a broad range of temperatures the only missing parameter for the temperature-dependent axion mass is the f a constant. For this several cosmological constraints exist [96]. To derive these model dependent assumptions for the axion couplings have to be made, e.g. using the DFSZ [97] or the KSVZ [98,99] model. The strongest lower bound is coming from the SN 1987A supernova event, suggesting f a 10 −9 GeV [96]. An upper bound is set by the known abundance of the cold dark matter. Assuming a random alignment for the axion field in the expanding Early Universe the axion zero mode obeys the damped oscillator equation It is initially completely dominated by the Hubble friction term. Only as the temperature drops and, thus, the topological susceptibility ∼ m 2 a rises such that m a ≈ 3H will the oscillatory dynamics start. Following the cooling rate of the universe from Eq. (4) and assuming an adiabatic evolution of the oscillations in the time-dependent axion potential the evolution of the axion density can be predicted. The energy density today mostly depends on the entropy ratio between today and the starting point CONF12 of the axion oscillations. Through this mechanism, the misalignment mechanism, the Universe is populated with axions, their density is strictly bound from above by the dark matter abundance [100].
Axions, if they exist, might provide only a small fraction of the dark matter, especially since the axionic string formation may also be non-negligible. Using the newly calculated topological susceptibility one can now calculate the cosmological axion bound. An axion mass of m a = 28(2)µeV represents a lower bound, in that case the axions from the misalignment mechanism could account for 100 % of dark matter, or, if the axion mass is m a = 50(4)µeV 50 % percent is accounted for. If one assumes that axions from the misalignment mechanism contribute between 1 and 50 % to dark matter, this range translates to a mass range of m a = 50 − 1500µeV [44].

Summary
In the past years the finite temperature lattice calculations have been extended to several new observables. Here we focused on the equation of state, and fluctuations and the topological susceptibility. The latter direction was the most demanding from the technical point of view, where severe discretization errors, an extremely low signal to noise ratio and long autocorrelation times had to be defeated. We hope that progress reported here will continue and soon the higher moments will be determined both in chemical potential and in the θ dependence.