Ideal topological gas in the high temperature phase of SU(3) gauge theory

We show that the nature of the topological fluctuations in $SU(3)$ gauge theory changes drastically at the finite temperature phase transition. Starting from temperatures right above the phase transition topological fluctuations come in well separated lumps of unit charge that form a non-interacting ideal gas. Our analysis is based on a novel method to count not only the net topological charge, but also separately the number of positively and negatively charged lumps in lattice configurations using the spectrum of the overlap Dirac operator. This enables us to determine the joint distribution of the number of positively and negatively charged topological objects, and we find this distribution to be consistent with that of an ideal gas of unit charged topological objects.

The presence of topologically nontrivial gauge field configurations is a peculiar feature of QCD that has important phenomenological consequences. Most recently this was highlighted by calculations to estimate the axion mass [1]- [3]. An essential ingredient of the calculation was the determination of the temperature-dependence of the topological susceptibility up to temperatures well above the QCD crossover temperature to the quark-gluon plasma.
At very high temperatures, fluctuations of the topological charge are strongly suppressed and occur in the form of localized lumps of action, carrying topological charge ±1. These objects are probably close in their shape and other properties to solutions of the classical Euclidean gauge field equations, i.e. instantons, or rather their finite temperature counterparts, calorons 1 [4]- [6]. Moreover, since at high temperatures, fluctuations of the topological charge are strongly suppressed, calorons (and anticalorons) are expected to form a dilute gas and their size is limited by the inverse temperature.
These properties of the caloron gas motivate the so called dilute instanton gas approximation (DIGA) which -in principle-makes it possible to calculate the temperature-dependence of the topological susceptibility perturbatively. However, the topological susceptibility determined in lattice simulations differs by an order of magnitude from the DIGA predictions even at temperatures as high as 5 − 10T c [1]- [3]. It is thus clear that at least one of the two assumptions that the DIGA is based on is not satisfied. These two assumptions, both of which are expected to be valid at high enough temperatures are: A1: The probability of one instanton occurring in a given volume can be calculated perturbatively in the semiclassical approximation. A2: The instantons gas is so dilute that interactions among (anti)instantons can be neglected, the gas of topological objects is an ideal gas.
Assumption A1 has been recently reconsidered, but despite the correction of the previously grossly underestimated uncertainty of the semiclassical one-instanton calculation, there is still at least a 3σ discrepancy between the lattice and the DIGA result for the topological susceptibility [7]. In the present paper we focus on assumption A2 and study interactions among topological objects in the quenched approximation of QCD, just above the finite temperature phase transition. Even apart from the axion problem, a full determination of instanton interactions just above T c is interesting in itself, as it can shed light on how typical gauge field configurations change as the system crosses into the high temperature phase. This might also help us better understand the chiral and deconfining transition in QCD with dynamical quarks.
To see how interactions among instantons could be detected, our starting point is a non-interacting instanton gas. In a free, non-interacting topological gas the number distribution of topological objects can be characterized with a single parameter, the topological susceptibility χ = Q 2 V , where Q = n i − n a is the topological charge (the number of instantons minus the number of antiinstantons), V is the space-time volume and . denotes the expectation with respect to the path integral. In an ideal topological gas all higher cumulants of the distribution can be exactly calculated in terms of χ. Any deviations from these ideal-gas cumulants are a result of interactions among topological objects.
In the recent literature several quenched lattice calculations of the lowest non-trivial cumulant 2 appeared [8,9]. The most precise calculation reports that even though above 1.15T c the value of B 2 is consistent with 1, its ideal-gas value, just above the phase transition, at 1.045T c it is still 1.27(7) [9]. For a more complete assessment of the situation, more information would be desirable about the distribution of the number of topological objects, beyond the first nontrivial cumulant of Q. However, higher cumulants of the distribution are notoriously hard to calculate, and even the full topological charge distribution can in principle miss subtle correlations among instantons and antiinstantons. Full information about that is contained only in the joint distribution of the number of instantons and antiinstantons.
The problem is that while there are well-established methods to compute the topological charge Q = n i − n a in lattice simulations, there is no easy way to determine n i and n a separately in lattice configurations 3 . In the present paper we introduce a novel method to compute n i and n a separately, and determine their joint distribution. Our method is based on the low-lying spectrum of the overlap Dirac operator. In particular, our main observation is that mixing instanton-antiinstanton zero modes constitute a distinct part of the Dirac spectrum close to zero, and can be reliably separated from the rest of the spectrum. Counting the number of these close to zero modes together with the exact zero modes of the overlap operator provides a way to determine not only the topological charge Q, but also the total number of topological objects n i + n a .
Besides yielding much more information than just the cumulant b 2 , our method has another advantage compared to previous studies. By construction it always gives integer numbers, and thereby avoids the ambiguities that plague the definitions of the topological charge based on gauge field operators. In that case, the values of the topological charge are not integers and need to be multiplicatively renormalized, and -if the full charge distribution is needed -also rounded to integers. Higher moments of the distribution, like b 2 are very sensitive to these ambiguities.
For the present study we used quenched lattice configurations generated at T = 1.045T c on lattices of temporal extension N t = 8 and aspect ratio 3 and 4. For both spatial volumes we determined n i and n a on 5k lattice configurations. In the smaller volume, we found that the number distribution of topological objects significantly deviated from the expectation based on a free non-interacting gas. In contrast, the larger volume showed no such deviation. Although in finite temperature lattice QCD an aspect ratio of 3 is usually considered safe in terms of finite (spatial) volume corrections, we show here that the unexpectedly large finite volume corrections are due to the proximity of the phase transition. We conclude that in quenched QCD already slightly above T c the number distribution of topological objects is consistent with that of a gas of free topological objects.
Let us first motivate the main tool used in our study, the separation of the bulk of the spectrum and the topology-related close to zero modes. It is known that in the presence of an isolated instanton (or antiinstanton) the Euclidean Dirac operator has an exact zero mode with chirality +1 (−1) [11]. In the field of a well separated instanton and antiinstanton the two would be zero eigenvalues split slightly and produce two complex conjugate eigenvalues. The splitting is controled by the spatial distance of the topological objects (relative to their size), as well as their orientation in group space. Generally the farther away the two objects are, the smaller the splitting is, and in the limit of infinite separation the splitting tends to zero [10]. In this way a dilute gas of topological objects is expected to produce not only |Q| = |n i − n a | exact zero modes, corresponding to the net topological charge, but also n i + n a − |Q| small Dirac eigenvalues, from the mixing of opposite chirality instanton and antiinstanton would be zero modes. Motivated by the instanton liquid model, we call the region in the spectrum containing these modes the Zero Mode Zone (ZMZ).
It should be already clear from the above discussion that as the temperature gets higher and topological fluctuations become sparser, the near zero modes of topological origin will be closer to the origin. At the same time, the low end of the bulk of the spectrum, the lowest non-topological modes, controlled by the Matsubara frequency, will move higher as the temperature increases. Therefore, at high enough temperature the ZMZ should be well separated from the bulk of the spectrum. In what follows, we will demonstrate that already slightly above the finite temperature phase transition the ZMZ can be reliably separated from the rest of the Dirac spectrum, provided a chirally symmetric Dirac operator, such as the overlap [12] is used.
Already in the early days of the overlap it was noticed that above T c , besides the expected exact zero modes, the overlap Dirac spectrum also contains an unexpectedly large number of very small close to zero eigenvalues [13]. This enhancement of the low end of the Dirac spectrum resulted in a spike in the spectral density, well separated from the bulk of the spectrum. This came as a surprise, since above T c the restoration of chiral symmetry would imply a vanishing spectral density at zero virtuality, due to the Banks-Casher relation. The spike in the spectral density was conjectured to contain mixing would be zero modes of instantons and antiinstantons. Subsequent work confirmed that this enhancement of the spectral density is neither a discretization, nor a quenched artifact [14]. More recently the appearance of this spike in the Dirac spectrum was speculated to signal a genuinely new state of strongly interacting matter, intermediate between the low temperature hadronic and the high temperature quark-gluon plasma state [16].
In the present work we analyze the statistical properties of the eigenvalues in this spike of the spectral density in a high statistics quenched SU (3) lattice study. We show that the statistics of these eigenvalues is to a high precision consistent with the assumption that they are produced by mixing instanton and antiinstanton would be zero modes. To this end we use quenched gauge ensembles generated with the Wilson gauge action at β = 6.09 and temporal lattice extension N t = 8. This corresponds to a temperature of T = 1.045T c , just above the finite temperature transition that in the quenched SU (3) case is a first order phase transition. For the detailed statistical analysis we used two ensembles of gauge configurations with spatial extension L = 24 and 32, both containing 5000 configurations. In addition, to check finite volume effects in the spectral density and the Polyakov loop distribution, we also had an ensemble of 600 configurations on a larger spatial volume L 3 = 40 3 . The negative mass parameter of the overlap Wilson kernel was set to M = −1.3, and two steps of hex smearing [17] were performed on the gauge links before inserting them into the Wilson kernel. The statistical analysis we report here was performed on the overlap Dirac eigenvalues of smallest magnitude with |λ|/T c < 2.0 for all ensembles.
Since we want to compare the statistics of small Dirac eigenvalues with that of non-interacting topological objects, we first summarize our expectations in such an ideal gas of topological objects. In the absence of any interaction among them, both the number of instantons n i and that of antiinstantons n a are expected to follow independent and identical Poisson distributions with mean V χ/2 proportional to the volume, where V is the four-volume of the system and χ will turn out to be the topological susceptibility. The joint distribution of n i and n a can be used to compute all the relevant physical quantities of this free topological gas in terms of the single parameter χ. In particular the topological susceptibility is where expectations like n i are understood to be with re- spect to the respective Poisson distribution. The topological charge distribution for Q ≥ 0 is where I Q are the Bessel functions of imaginary argument. Due to time-reversal symmetry the distribution is symmetric, P (Q) = P (−Q). Another interesting quantity to consider is the distribution of the total number of topological objects n = n i + n a , which is simply a Poisson distribution with mean V χ.
Let us now confront the lattice data with these expectations. In Fig. 1 we show the spectral density of the overlap Dirac operator on the previously mentioned lattice ensembles. Although in the spontaneously broken phase of the pure gauge theory that we consider here, the three Polyakov loop sectors are equivalent, the spectrum of the Dirac operator and the pattern of chiral symmetry breaking/restoration is known to be strongly dependent on the Polyakov-loop sector [6], [18]- [21]. Since in the present work we want to study features that are expected to be at least qualitatively present in QCD with dynamical quarks, here we restrict the analysis to configurations in the physical ReP > 0 Polyakov loop sector. This is the only one that would be present if dynamical fermions were to be included. In fact, even the slightest explicit breaking of the Z(3) symmetry by dynamical fermions with arbitrarily large, but finite mass would force the system into the real Polyakov-loop sector. The enhancement of the spectral density near zero is clearly seen in Fig. 1. We note that the exact zero eigenvalues are not shown here, they would appear as a delta function exactly at zero.
Counting the number of zero eigenvalues allows us to compute the topological susceptibility χ, as well as the distribution of the topological charge. In Fig. 2 we compare the distribution obtained in the lattice simulation with the one expected in a free topological gas with susceptibility χ. This is essentially a one-parameter fit of the function in Eq. (4), the fit parameter being V χ and the chi squared per degree of freedom of the fit turns out to be 0.85.
Encouraged by the good agreement between the lattice data and the free topological gas, we assume that the exact zero modes and the small Dirac eigenvalues, up to a point λ ZMZ in the spectrum, are the eigenvalues associated to the topological objects. In this way, by counting them we count the number of topological objects n present in the gauge field. To make this picture consistent, we have to choose λ ZMZ such that n = V χ, as predicted by Eq. (5) for an ideal topological gas. Requiring this, results in λ ZMZ a = 0.045(6) 4 , which turns out to be in the depleted region of the spectral density, separating the spike at zero from the bulk (see Fig. 1). This shows that the ZMZ is indeed well separated from the bulk of the spectrum.
We can now identify the total number of eigenvalues in the zero mode zone, the ones with |λ| < λ ZMZ (including the exact zero modes) with n, the number of topological objects. Counting the eigenvalues in the ZMZ configuration by configuration, we obtain the distribution of n and in Fig. 3 we compare it with the one expected in a gas of noninteracting topological objects, given by Eq. (5).
We emphasize that at this point no fitting is involved, since the only parameter of this distribution, V χ, had already been determined independently from the charge distribution. We do not find a significant deviation from the free topological gas distribution, as the chi squared per degree of freedom of the deviation is 0.62. We emphasize that the fact that the distribution of the number of eigenmodes of magnitude smaller than λ ZMZ exactly reproduces the distribution we expect from eigenmodes linked to free topological objects is rather nontrivial. By choosing λ ZMZ in the above described manner, we only fixed the expectation of the distribution, and it follows the expected analytical form over three orders of magnitude, the whole range where we have data. This shows not only that the topological objects are indeed non-interacting, but also that already at this temperature the zero mode zone can be reliably separated from the bulk.
So far we determined λ ZMZ from the requirement that the total number of eigenvalues below λ ZMZ (including the zero modes) be consistent with the topological susceptibility obtained by counting the zero modes only. But is λ ZMZ really a special point in the spectrum? To answer this question we chose different cuts in the spectrum and checked how close the distribution of the number of eigenvalues below these cuts are to a Poisson distributions. To this end we determined the distribution of the number of eigenvalues as a function of the cut and for each value of the cut plot the chi squared per degree of freedom of the deviation of the best fit Poisson distribution from the data. In Fig. 4 we show the results. It can be clearly seen that there is a range of cuts 0.3 < λ cut /T c < 0.6, where the distribution is compatible with Poisson, and the previously independently determined λ ZMZ happens to be in the middle of this range. This further supports that the valley of the spectral density containing λ ZMZ indeed separates the topological modes from the bulk of the spectrum. This valley, however, is rather wide and even though the spectrum there is quite depleted and not many eigenvalues are contained there, the question arises as to how sharply λ ZMZ is defined within the valley. Given the present data set, this question cannot be unambiguously answered. It is possible that with much larger statistics the chi squared test presented in Fig. 4 would further limit the acceptable range for λ ZMZ . Another possibility is that in the continuum limit the spectrum could become more depleted in the valley, even a gap could appear there. In that case the exact location of λ ZMZ within the gap would not be important. Finally, it is also possible that even in the continuum limit and with arbitrarily large statistics there would still be some small ambiguity in separating the topological modes and the bulk. To explore these possibilities further would require more extensive simulations that are out of the scope of the present work. The finite temperature SU (3) transition is a first order phase transition, so the correlation length does not diverge, however, large finite volume corrections cannot be excluded. To assess their importance, we repeated the analysis in a smaller volume with linear extension L = 24. In that case we found significant deviations from the expected free topological gas behavior. The resulting chi squared per degree of freedom was 1.99 and 6.29 in the case of the charge distribution and the distribution of the total number of topological objects, respectively.
To understand finite volume corrections in the vicinity of a phase transition, it is instructive to look at the volume dependence of the distribution of the order parameter. In Fig. 5 we show the probability distribution of the magnitude of the Polyakov loop, the order pa- rameter of the quenched finite temperature transition.
Besides the widening of the distribution, expected for smaller volumes, the data for L = 24, 32 lattices also show an unusual enhancement of smaller values of the Polyakov loop. The reason for this is that in the high temperature phase the Z(3) center symmetry is spontaneously broken and the system randomly chooses one of the three Z(3) sectors. However, in a finite volume tunneling among the sectors is still possible, the tunneling probability is enhanced in smaller volumes, and configurations in the process of tunneling have small magnitudes of the Polyakov loop. As also seen in Fig. 5, in larger volumes these tunneling states get suppressed, however, in smaller volumes they can still give significant contributions to physical quantities, resulting in large finite-size effects. To see, how these tunneling states can affect the topological charge, we looked at the correlation between topology and the Polyakov loop. The simplest quantity to study is the topological susceptibility. We computed its dependence on the Polyakov loop by restricting the averaging of Q 2 to configurations with Polyakov loop magnitudes in intervals of length 0.01. The results for the L = 24 and 32 ensembles, shown in Fig. 6, reveal a strong dependence of the susceptibility on the Polyakov loop. The previously seen enhanced contribution of the tunneling region (small Polyakov loop), where the susceptibility is larger, will thus result in significantly larger topological susceptibilities in smaller volumes. To have a feeling about the relative importance of the enhanced region, we note that the probability that |P | < 0.08 is 0.11, 0.028 and 0.004 for the lattices with linear spatial size size L = 24, 32 and 40 respectively. For the two ensembles shown in Fig.  6 we also indicated the overall susceptibilities and their uncertainties with the horizontal stripes. Since on the L = 24 lattices even the susceptibility suffers large finitesize effects, it is not surprising that the same occurs for the distribution of the topological charge and the number of topological objects.
We would also like to comment on the apparent discrepancy between our results and those of Ref. [9] who found a significant (27(7)%) deviation of the B 2 coefficient from its value (1.0) expected in a non-interacting instanton gas. In fact, our data yields B 2 = 1.35(41), which is compatible with that of the above reference. However, judging from their much smaller uncertainty, their statistics could be more than an order of magnitude larger than ours. Since it is based on the overlap spectrum, our method is computationally much more expensive, and in the present study we could not compete in the statistics, but our method has two advantages. First, it necessarily yields integer charges and avoids the large ambiguity in B 2 due to any small random fluctuations of the topological charge around integers and the possibly necessary normalization and rounding of the charges. Secondly, our method allows for a full determination of the joint distribution of instantons and antiinstantons. It would be interesting to repeat our calculation using a much larger statistics and to see if there is any deviation in this distribution from that of the number of noninteracting instantons.
In the present paper we used a novel way to compute the joint distribution of the number of topological objects in lattice simulations. We showed that right above the critical temperature of pure SU (3) gauge theory the distribution is consistent with the one expected in an ideal gas of non-interacting charges of unit magnitude. It is remarkable that while below the phase transition topological fluctuations form a dense medium without easily identifiable individual lumps [22], right above the phase transition an ideal gas of well separated topological lumps emerges. Our result also implies that the most likely explanation of the large discrepancy between the lattice and DIGA based calculation of the topological susceptibility is that the topological lumps we found are not close enough in shape to ideal calorons to warrant a semiclassical treatment.
We expect that -at least on a qualitative level-this picture of the topological fluctuations that we found in the quenched case carries over to QCD with dynamical quarks. However, the fermion determinant might introduce some interaction even among well separated topological lumps, but to study that one would need to use a chiral Dirac operator also for the simulation of the sea quarks.