The topological properties of QCD at high temperature: problems and perspectives

Lattice computations are the only first principle method capable of quantitatively assessing the topological properties of QCD at high temperature, however the numerical determination of the topological properties of QCD, especially in the high temperature phase, is a notoriously difficult problem. We will discuss the difficulties encountered in such a computation and some strategies that have been proposed to avoid (or at least to alleviate) them.


Introduction
Since the discovery of the instanton solution [1] and the introduction of the θ-vacuum [2,3] there has been constant interest in the study of the intrinsically non-perturbative aspects of gauge theories related to the θ-dependence. Such an interest was motivated by phenomenological and theoretical reasons: from the theoretical point of view θ-dependence is a new knob that can be used to investigate the properties of strongly interacting gauge theories, from the phenomenological side it is connected with the anomalous U A (1) symmetry of massless QCD and the CP invariance of strong interactions.
The experimental fact that strong interactions are invariant under the discrete CP symmetry is indeed related to the unnaturally small value of the (dimensionless, renormalization group invariant) θ parameter entering the QCD lagrangian, with the best present estimate |θ| 10 −10 being obtained from measures of the neutron electric dipole moment [4]. One of the most appealing way of explaining this small value of θ is to assume it to be in fact exactly zero due to some new symmetry. This approach was initiated by Peccei and Quinn [5,6] and, soon after its proposal, this mechanism was shown to imply the existence of a new light pseudoscalar particle [7,8], the axion, whose properties are strictly related to the behaviour of the effective potential for small θ values.
More or less during the same years the theoretical study of θ dependence in QCD was systematized and, in particular, it was realized that two different behaviours are to be expected in the low and in the high temperature regimes. At T = 0 chiral perturbation theory can be used to compute the effective potential as a function of the θ value [9,10], obtaining at leading order for the physical case of two light flavours the result In the very high temperature regime one can instead use perturbation theory around instantons to show that the instanton density gets lower and lower as the temperature increases, and this implies for the free energy density the functional form [11] f (θ, where the subscript "DIGA" stands for Dilute Instanton Gas Approximation. When N f light flavours are present, the coefficient χ DIGA (T ) scales with the temperature as [11] χ DIGA (T ) ∝ m N f T 4− 11 Finally, in [12][13][14], the axion field was recognized to be a good cold dark matter candidate and it was realized that this implies an overclosure bound: using the knowledge of the θ dependence of QCD at high temperature we can estimate the number of axions produced in the ealry universe and thus the axion density in today universe. From the fact that this density has to be smaller than the dark matter density we obtain constraints on the axion interactions and, in particular, a lower bound for its mass (modulo some cosmological hypotheses), see [15] for more details. In order to have reliable bounds (or estimates, if we assume all the dark matter to be composed of axions) we need to know with reasonable accuracy the θ dependence of QCD at finite temperature and, in particular, across the deconfinement/chiral restoration transition. The same information is also needed if we are interested in studying the interplay between topology and color confinement in gauge theories.
Since solid analytical results for the θ dependence are available only for the cases of very low temperature (from chiral effective field theory) and asymptotically high temperatures (from perturbation theory), Lattice QCD is the only approach available to systematically investigate the phenomena taking place for temperatures of the order of the critical temperature T c . Some of the points we would like to clarify are: • up to which temperature is it reasonable to use the predictions of chiral perturbation theory?
• what happens for T close to T c ?
• for which temperature we can start using DIGA results?
In order to answer these questions it is convenient to use a general parametrization of the θ dependence of the free energy density f . Assuming f (θ, T ) to be analytic in θ and using the simple property f (θ, T ) = f (−θ, T ) we can write the free energy density in the general form [16] Since θ enters the Euclidean lagrangian only through the coefficients χ(T ) and b 2n (T ) defined by Eq. (4) can be related to the cumulants of the topological charge Q computed at vanishing θ; it is for example easy to show that the topological susceptibility χ(T ) and the coefficient b 2 (T ) are given by the expressions where V is the four dimensional volume of the lattice in physical units. At T = 0 we have (for physical up and down quark masses) the next to leading order χPT values χ(T = 0) 1/4 = 75.5(5) MeV and b 2 (T = 0) = −0.029(2) (see [10]), while in the DIGA regime it is immediate to get from Eq. (2) the values b 2 = −1/12 and b 4 = 1/360.
To investigate the θ dependence of QCD by means of lattice simulations we have to evaluate, using Eq. (6) or equivalent forms, the coefficients χ(T ) and b 2n (T ) appearing in the general parameterization of the free energy density in Eq. (4). As usual, the pure Yang-Mills theory is a convenient test bed to inquire in a computationally simpler setup the level of accuracy that one can hope to reach in QCD. The numerical study of the topological properties of Yang-Mills theory at finite T has quite a long story and has reached solid results on basically all the open fronts. It is indeed known that • the topological susceptibility is practically independent of the temperature up to deconfinement (see e.g. [17] for an early study), • close to T c the functional form of the θ dependence abruptly changes and for temperatures slightly above deconfinement Eq. (2) well describes data (see [18] and later improvements in [19,20]), • the DIGA-like scaling of the topological susceptibility χ(T ) ∝ T −7 starts for the S U(3) case at surprisingly low temperatures [21][22][23] but DIGA predictions underestimate χ(T ) by about an order of magnitude in the temperature range T c T 4T c .
It is important to note that these results are not easy to obtain even for the pure glue case. The obstacles encountered are similar to the ones that are found also in the QCD case, however in Yang-Mills theory, where more efficient algorithms and higher statistics are available, it is possible to proceed by brute force. The two main problems encountered are the so called freezing of the topological charge and the presence of particularly severe finite size effects at high temperature.
What goes under the name of "freezing of the topological charge" is the fact that the autocorrelation time of the topological charge Q quickly becomes very large as the continuum limit is approached, until ergodicity is lost and reliable simulations cannot be performed any more. Since Q is an integer (or almost integer, depending on the discretization adopted) number, in simulations performed at (too) Traj. Traj. fine lattice spacing the topological charge stays constant during all the simulation time, i.e. it is frozen. An example of this behaviour is shown in Fig. 1(left). The fact that the autocorrelation time of Q gets larger and larger as the continuum limit is approached is quite easy to understand: since Q is a topological invariant, to change its value it is necessary to pass through configurations that are practically discontinuous and these configurations are more and more suppressed as we get closer to the continuum limit. Given this simple argument we could guess the autocorrelation time of Q to scale approximately as the exponential of the correlation length in lattice units. This scaling form is compatible with numerical results obtained in QCD [24,25] and in other systems with topologically stable classical solutions [26][27][28][29][30], however it is notoriously difficult to distinguish numerically an exponential behaviour from a power law behaviour with very large exponent. Despite the uncertainties related to the precise functional form of the critical slowing down of the topological charge, it is clear that it is much more severe than the slowing down observed for "non topological" observables, and this puts strong constraints on the lattice spacing values that can be used in simulations.
In the high temperature phase the freezing problem is still present (although it appears at different values of the lattice spacing for different temperatures) however in this case there is also a new source of trouble: since the topological susceptibility goes to zero at high temperature, by keeping the physical volume fixed and increasing the temperature we are reducing the typical amount of topological charge that is present in the lattice, which is given by χ(T )V (see Eq. (6)). When this number becomes much smaller than 1 we are effectively sampling only the Q = 0 topological sector of the theory and it is very difficult to reliably estimate both χ and the b 2n coefficients. In Fig. 1(right) an example of this phenomenon is shown in a case in which it is still not extreme, i.e. in which also the Q = ±1 sectors are occasionally explored.
From the practical point of view this problem is analogous to the freezing problem, however these two problems have deeply different origins: freezing is a purely algorithmic issue and it is at least conceivable that an algorithm exists that completely solves it 1 , the other problem is instead of physical origin and it is related to the fact that at high temperature, on a finite volume, the statistical weight of the Q = 0 configurations becomes overwhelming in the path-integral. As a consequence freezing can happen for every value of the topological charge, while in the high temperature phase only the probability of the Q = 0 sector is enhanced.
The two problems mentioned so far are common to both QCD and pure gauge theories, however when dynamical light fermions are present a new problem emerges: the convergence to the continuum limit of the topological susceptibility is extremely slow and very fine lattice spacings are needed to extract the continuum limit, see [31][32][33]. The origin of this slow convergence rate to the continuum is the strong connection between the θ dependence and the chiral properties of the dynamical fermions, that can be seen as a consequence of the index theorem and of the chiral anomaly. Since in most lattice discretizations the chiral symmetry is at least partially broken, large lattice artefacts appear in the θ dependence and, in particular, in the topological susceptibility.
To summarize the computational difficulties: to obtain reliable results for the θ dependence we have to perform simulations at very fine lattice spacings, since otherwise huge lattice artefacts are present. On the other hand we cannot use very fine lattice spacings since, beyond the usual computational burden, we also have the topological freezing problem. All this remains true in the high temperature phase, where a strong suppression of the Q 0 sectors is also present. Just to make things even worst the topological susceptibility goes to zero as we lower the quark masses, and thus the high T finite size effect worsen when using dynamical light fermions.
Of the three problems discussed in the previous section the one that attracted most of the activity in the past was the freezing problem, likely because it appears not only in QCD but in all field theories (or even quantum mechanical models) whose configuration space can be decomposed, in the continuum, into homotopy classes. This property makes the freezing problem amenable of being studied using models of increasing computational complexity before attacking the QCD case; typical examples of systems used in test studies are: quantum mechanics on a circumference, two dimensional CP N models and Yang-Mills theories. In the years several methods have been developed to cope with the freezing problem and two main approaches can be found in the literature: the goal of the first group of algorithms is to reduce the autocorrelation time of the topological charge, while methods of the second group aim at extracting as much information as possible from frozen simulations.
One of the most widely used method of the first group is the use of open boundary conditions [34]. The idea of the method is quite intuitive and it is related to the fact that with periodic boundary conditions (i.e. effectively no boundaries) the only possible way of changing the topology is to pass through discontinuous configurations having very small probability in the path-integral; if instead we use open boundary conditions there can be an inflow/outflow of topology from the boundary and the autocorrelation time of Q is expected to scale as ∝ a −2 [35]. Simulations adopting open boundary conditions has been by now used mainly to study the T = 0 case, in which it is natural to "open" the torus in the temporal direction, however it is in principle possible to use the same idea also at finite temperature, using open boundary conditions in a spatial direction [36]. An obvious drawback of this method is however the fact that we loose translation invariance 2 and we have to study only the part of the lattice that is far enough from the boundaries not to be significantly affected by them. In particular, expressions like Eq. (6) cannot be used anymore, and the topological susceptibility has to be computed as the integral of the two point function of the topological charge density q(x), with the two point function being estimated in the bulk of the lattice. The b 2n coefficients can be in principle computed in an analogous way, using the 2n−point functions of q(x), however this becomes increasingly difficult for larger n values.
Another method to improve the decorrelation of the topological charge is the use of extended space approaches, like e.g. tempering methods. These have been introduced in the study of spin glass systems and the idea is to mix together in a stochastically exact combination simulations performed at different values of the parameters, in such a way that the slowly decorrelating runs are speed up by the quickly decorrelating ones. Two different practical implementations of this general philosophy are simulated tempering [38] and parallel tempering [39]. In the spin glass case the value of the autocorrelation time increases at low temperature and the parameter that is used in the tempering procedure is thus temperature; for the case of the freezing problem a natural choice is the lattice spacing. This approach has been tested in toy models with promising results (see [40] for CP N models and [30] for a quantum mechanical model), but a systematical investigation for the case of lattice gauge theories has not yet been performed. A different application of the tempering approach has been advocated in [41], where it has been shown (for the case of CP N models) that parallel tempering can be efficiently used to connect simulations performed using, instead of different parameter values, different boundary conditions, in order to reduce the finite size effects of the open boundary conditions while retaining some of their advantages.
A third class of approaches that can be used to reduce the autocorrelation time of Q in simulations is to sample, instead of the original distribution, a different one that is computationally less problematic and finally reweight the results; this is the idea of the multicanonical algorithm, introduced in [42] for the study of strong first order transitions. Multicanonical simulations are usually carried out in such a way that configurations are sampled with a probability proportional to the inverse of the density of states (of the original system), in order to obtain a flat histogram in the modified system. As it always happens, there are several ways of implementing this general idea. The original proposal [42] was to use a rough estimate of the density of states, since the algorithm is stochastically exact anyway, however it was later realized that some improvements are possible on this basic approach. A possibility that can be thought as a dynamical implementation of the multicanonical idea is metadynamics [43], in which an history dependent biasing potential is added to the system and updated in such a way that, as the simulation goes on, the biasing potential converges to the logarithm of the density of states. This idea has been successfully applied to the case of CP N models in [44] and we will discuss some more details of metadinamics in the next section, where the application to the QCD case will be analyzed. A more radical approach is to precisely estimate the density of states [45,46] and then compute observables without the need of a Monte Carlo. First results for the topological susceptibility obtained in this way have been presented in [47].
As anticipated before, there are also algorithmic developments that go in the direction of obtaining as much information as possible from simulations in which the topology is frozen. An idea that is around since some time [48] but that recently draw new attention is the use of sub-volume measurements: even if the global topological charge is fixed we can try to extract information from its local fluctuations. In [49] a Maximum-Likelihood method was advocated to estimate the parameters of the θ dependence, while the approach proposed in [50] is more direct but, in its original formulation, it is limited to the case in which the b 2n coefficients vanish. It is nevertheless reasonable to expect that, by sacrificing part of its simplicity, a generalization of this algorithm to the case b 2n 0 can be found. Another interesting possibility, that was originally introduced for the case of overlap fermions simulations, is to extract the topological susceptibility from finite size effects of simulations performed at fixed Q [51,52]. An algorithm that uses also simulations at fixed topology is the one introduced in [33,53] (see also [54]): if for some values of the parameters (e.g. at coarse lattice spacing) it is possible to estimate the relative weights of the different topological sectors, then it is possible to propagate this information to other parameter values using an analogous of the thermodynamic integration for Z Q /Z 0 , where Z Q is the partition function at fixed Q value. From this information is then possible to evaluate both the topological susceptibility and the b 2n coefficients, the only possible difficulty being the fact that a large number of simulations could be needed if Q 2 1. As previously explained in Sec. 2, in simulations performed in QCD with dynamical light flavours the freezing of the topological charge is a particularly severe problem because of the fact that the topological susceptibility converges to its continuum limit very slowly. An approach that is complementary to the development of algorithms to cope with the freezing problem is thus the study of techniques to speed up the approach to the continuum limit for topological observables. Once realized that the slow convergence to the continuum can be related to the lattice breaking of the chiral symmetry, it is natural to investigate the possibility of correcting a posteriori this breaking in the numerical results to improve the convergence.
In the T = 0 case we can rely on chiral perturbation theory at finite lattice spacing in order to have a better understanding of the observed numerical behaviour. Staggered fermions have been up to now the most widely used discretization in studies of the θ dependence 3 and a staggered chiral perturbation analysis of the topological susceptibility has been carried out in [57]. The outcome of this analysis is that, at non-vanishing lattice spacing, the dominant violation of the expected chiral behaviour of the topological susceptibility can be described by using the taste-singlet pion mass m π,ts instead of the physical pion mass m π . The idea of correcting the numerical results by a factor related to the taste splitting was introduced in [31] and more thoroughly investigated in [33] and it was shown to be very effective in improving the convergence to the continuum limit of the topological susceptibility at zero temperature.
The finite temperature case is more problematic, since there are no solid theoretical results available for generic temperatures. In [31] it was proposed to use the dimensionless ratio χ(T, a)/χ(T = 0, a), with the idea that a large part of the lattice artefacts would cancel in this ratio. This appeared to be the case for the lattice spacings explored in [31], however it was later shown in [32,33] that the problem of the slow convergence to the continuum is not solved using this recipe: when going to still smaller lattice spacings the ratio χ(T, a)/χ(T = 0, a) starts again to show a strong dependence on the value of a. In the Dilute Instanton Gas Approximation each instanton (or anti-instanton) is associated to (at least) a zero mode of the Dirac operator and, since instantons are weakly interacting by hypothesis, the mass dependence of the partition function in the chiral limit is dominated by these isolated zero modes. The approach proposed in [33] (see also [54]) consists in reweighting gauge configurations in order to fix this behaviour. Let's take the example of a Q = 1 configuration: there should be at least a zero mode of the Dirac operator 4 however, due to the explicit lattice breaking of the chiral symmetry, only a quasi zero mode with eigenvalue ∆ will be present, in which case a weight m m+∆ N f is assigned to the configuration, where m is the light quark mass. Using this approach the value of the topological susceptibility that is obtained in the high temperature phase is reduced by more than an order of magnitude and the dependence of the result on the lattice spacing is very mild.
Let us close this section with some comments on the finite volume problem at high temperature. Since for large enough temperature the system is practically always in the charge zero sector, all the methods that have been developed to extract information from frozen simulations are in principle applicable also to this case. In particular the "thermodynamic integration" method of [33,53] has been originally introduced in this context and not as a solution to the freezing problem. Also the approaches that have been broadly classified before as "multicanonical like" (i.e. multicanonical, metadynamics and density of states) can be profitably used at high temperature, see the next section and [30].

Some results from metadynamics
In this section we will present the outcome of some tests carried out in collaboration with M. D'Elia, G. Martinelli, F. Negro and F. Sanfilippo to investigate the possibility of using metadynamics in the high temperature phase of QCD.
Metadynamics has been introduced to efficiently sample the configuration space of models with several free energy minima for which we have been able to identify (at least approximately) the slow degrees of freedom. From this point of view metadynamics is less versatile than e. g. tempering methods, however, when it can be applied, it is expected to be more efficient, since it explicitly uses more information on the system. The general philosophy of metadynamics is the following [43] (see [58] for a review) 1. identify the "slow variable" φ, typically called "collective variable" in the original literature (note that several slow variables can be used together) 2. introduce in the action an external history-dependent potential V(φ, τ) (called metapotential) for the slow variable, where τ is the Monte Carlo time 3. during the update modify V(φ, τ) in such a way that the distribution of φ is (on average) constant in the interval of interest [φ min , φ max ] 4. at the end of the simulation reweight data using V(φ, τ) , that is an estimator of the free energy at fixed φ.
The point 3) is typically implemented by increasing the value of the metapotential corresponding to the value of φ that is currently in use, in order to escape from free energy minima and explore uniformly the interval [φ min , φ max ]. As noted in the previous section this approach can be intuitively thought as a dynamical multicanonical update, since we start from a very rough estimate of the metapotential (like e.g. a flat distribution) that is then self-consistently improved during the simulation. Metadynamics, being a random walk biased toward exploring new values, is expected to converge faster than the multicanonical update, that corresponds to a simple random walk; on the other hand the algorithm does not use anymore a stationary Markov chain and, as a consequence, it is non trivial to prove its stochastic exactness [59].
In order to apply metadynamics we have first of all to identify the appropriate collective variable to be used. A natural choice would be the topological charge Q, however this would not work. If we use a discretization that gives an exactly integer value for Q, like the geometric or fermionic definitions, the metapotential induces no force in the evolution unless there is a topology change in the trajectory; as a consequence the "bias" would not be effective in the dynamics. Also the use of a definition of the topological charge that is not integer by construction, like the one obtained by using the clover discretization after having applied cooling or gradient flow, would not be appropriate. Indeed this topological charge is no more exactly integer, but it is still so strongly peaked at integer values that it would not sufficiently bias the dynamics. A pictorial representation of this fact is reported in Fig. 2. In Fig. 2(left) a free energy profile is sketched for which the application of metadynamics is expected to be useful: updating the metapotential we will gradually fill the local minima and the whole configuration space will then be sampled. If we use the topological charge Q as collective variable the situation is however more similar to Fig. 2(right): even using metadynamics it would be difficult to fill the minima and escape the freezing.
From the previous analysis it emerges that a good collective variable should be something that has good correlation with the topological charge Q but that can assume also far from integer values. Our choice was to use what could be called an under-smoothed topological charge, i.e. the result obtained by computing the clover discretization of the topological charge after just a few smearing steps. We will denote this variable by Q and, more precisely, it is defined as the value obtained by using the clover discretization of the topological charge after having applied 20 stout smearing steps [60] to the configuration, with smearing parameter α = 0.1. While a priori other smoothing procedures can  be adopted as well in the definition of Q , the use of a differentiable smearing is fundamental for the implementation of metadynamics in a standard HMC code.
A test run was performed by applying metadynamics to the collective variable Q in a simulation with N f = 2 + 1 flavours of stout smeared staggered quarks with physical quark masses and Symanzik tree level improved gauge action. The physical parameters adopted were the same as for the time history shown in Fig. 1(right), corresponding to a 0.40 fm and T 310 MeV, and metadynamics was used to improve the sampling of the interval [Q min , Q max ] = [−3 : 3], i.e. the metapotential takes the constant values V(Q = 3, τ) for Q > 3 and V(Q = −3, τ) for Q < −3.
In Fig. 3 the time histories of Q and Q in the metadynamics run are shown, together with their histograms, and several things can be noted: first of all it is clear that there is a good correlation between Q and the topological charge Q measured after 70 cooling steps 5 , moreover the histograms show that Q is strongly peaked on nearly integer values, while this is not the case for Q , whose histogram is nearly flat due to metadynamics. These were the requirements for the biasing potential to be effective and indeed we can see that, notwithstanding the temperature of about 310 MeV and the small lattice spacing (a 0.04 fm), all the topological sectors in the range [−3, 3] are explored and no signal of freezing is present. The overhead introduced by metadynamics is not negligible, since the time required to evolve a trajectory when metadynamics is used is about three times the one required in a standard HMC simulation; given the obtained results this seems however totally worth it.

Conclusions
The study of the topological properties of non abelian gauge theories, and QCD in particular, using lattice simulations has by now a long tradition, dating back to early eighties [61]. At low temperature chiral perturbation theory can be used, for very high T perturbation theory and instanton calculus are available, but the lattice approach is still the only tool that can be used to investigate, from first principles and in a systematically improvable way, the θ dependence of QCD for generic values of the temperature.
The fact that lattice simulations can (or have to) be used does not mean that the numerical determination by lattice methods of the θ dependence is an easy task, and in fact it is a surprisingly hard one. Mainly for this reason most studies focused in the past on the pure glue case, which is also problematic but for which a brute force approach can be pursued with success. The recent phenomenological interest in axion properties triggered new lattice studies of θ dependence in QCD with physical quark masses, that came together with the developments of new algorithms.
The main problems that are encountered when trying to estimate the θ dependence of QCD with physical quark masses are the exceptionally large autocorrelation time of the topological charge at small lattice spacing (the so called "freezing problem"), the fact that at high temperature the weight in the path-integral of the Q = 0 sector becomes overwhelming, the exceedingly slow convergence to the continuum limit of the topological susceptibility (and maybe also of higher cumulants of the topological charge).
The freezing problem was the one that attracted most of the algorithmic developments, likely because it is not only typical of QCD but it appears in all theories with topologically stable classical solutions. Nevertheless, in the last couple of years, several new algorithms have been developed to cope with all the previously mentioned problems and we are now in position to crosscheck the different approaches that have been developed, in order to systematically inquire the possible presence of systematics and provide solid results for the QCD case.