Towards a semi-classical description of QCD vacuum around $T_c$

We study the vacuum topology of 2+1 flavor QCD above the chiral crossover transition, at $T \lesssim 1.2~T_c$, on lattices of size $32^3\times 8$. Since overlap fermions have exact chiral symmetry and an index theorem even on a finite lattice, we use them to detect the topological content of gauge fields generated using domain wall fermion discretization for quarks. We further use different periodicity phases along the temporal direction for the valence overlap quarks, which allows us to probe different topological structures present in the gauge field ensembles, through its zero modes. This procedure provides strong evidences that fermion zero-modes can be quantitatively understood to arise due to different species of instanton-dyons. We estimate their relative abundances from the Dirac-eigenvalue density and resolve the so called"topological clusters"via multi-parameter fits to their density, providing therefore an understanding of the interactions between instanton-dyons. The typical separation between dyons we obtain, is $\sim 0.3$ fm. Surprisingly, it emerges out from this study that a semi-classical description of the fermionic zero modes in the QCD vacuum is quite accurate just above $T_c$.


I. INTRODUCTION
With the increasing sophistication of non-perturbative lattice gauge theory techniques, we now know that the (nearly exact) chiral symmetry in 2+1 flavor Quantum Chromodynamics (QCD) is restored via a continuous and smooth crossover transition [1][2][3]. The chiral crossover transition temperature was recently quantified to an unprecedented accuracy [4] and is T c = 156.5 ± 1.5 MeV. The observables sensitive to chiral symmetry breaking and confinement show consistent behavior in the crossover transition, implying a simultaneous occurrence of these two phenomena in QCD with quarks in the fundamental representation of the gauge group [5]. The question that still remain unresolved is the precise topological identity of the microscopic degrees of freedom responsible for these two physically distinct phenomena in QCD.
With the motivation of explaining confinement in QCD, Polyakov introduced instantons. In the 2+1 D case he showed in his seminal work [6] that a gas of instantons (monopoles in this setting) produces the famous confining potential. However in the (most important) 3+1D case, the ensemble of instantons were not confining due to the fact that their interaction is not long-range. * Electronic address: rlarsen@bnl.gov † Electronic address: sayantans@imsc.res.in ‡ Electronic address: edward.shuryak@stonybrook.edu The instanton liquid model (ILM) [7] used the instanton density and the average radius of the instantons as inputs. With their values known, one can show that the instanton-induced t'Hooft Lagrangian is strong enough to break SU A (N f ) chiral symmetry, like it was the case for the Nambu-Jona-Lasinio (NJL) model [8]. Further studies of interacting instanton ensembles predicted different hadronic correlation functions within ILM, which are in fairly good agreement with the results from QCD inspired models and lattice [9]. This model explained not only the SU A (N f ) chiral symmetry breaking and pions, but also large violations of the U A (1) symmetry, e.g. the large mass of the η meson. However confinement could not be explained within this framework. The (renormalized) Polyakov loop is traditionally used as the order parameter of the confinement transition in pure gauge theory. In order to understand both deconfinement and chiral symmetry breaking, one needs to find topological solitons sensitive to the eigenvalues of the Polyakov loop. At finite temperatures, such solitons with a non-trivial holonomy have been shown to exist in [10][11][12][13]: they are called instanton-dyons or instantonmonopoles. For simplicity, we henceforth refer to them as "dyons". Thermodynamic ensembles of dyons were studied analytically [14][15][16][17] via mean-field approximation, as well as numerically [18][19][20]. These works could explain (near simultaneous) occurrence of both deconfinement and chiral phase transitions in QCD-like theories. Unlike instantons, the instanton-dyons carry not only topological but also electric charges, thus back-interacting with the Polyakov loop inducing confinement [21,22].
They also carry magnetic charges and therefore Dirac strings, which explains why their topological charges are not quantized as integers.
Whether this is indeed the right explanation or not, can be directly tested using lattice gauge-field configurations. One way to do so is by using over-improved cooling technique [23,24]. The other one, which we use, is via studying the properties of the fermionic zero and nearzero modes, with modified temporal periodicity phases as a tool [25]. Studies along these lines, have been performed by Mueller-Preussker, Ilgenfritz and collaborators, see for e.g. [26,27] which provided tantalizing evidences for the presence of dyons also in finite temperature QCD.
The aim of this work is to provide quantitative description of these Dirac eigenstates in terms of superposition of the dyons. Using sophisticated lattice techniques we aim to isolate and identify different species of dyons and understand the nature of interactions between them. For the first time, fermions with exact chiral symmetry on the lattice, the so-called overlap fermions, are used for this aim. Our initial results have been published in Ref. [28]. Even though computationally expensive, such fermions are ideally suited to probe the topological content of the gauge configurations, which, in our case, are generated with the domain-wall fermion discretization. By varying the temporal periodicity phases of the (valence) fermions, we track all fermionic zero modes. We show that not only the index theorem works as expected, but the objects responsible for the zero modes are indeed nothing but different species of dyons. Further nontrivial tests of this picture is a comparison of the quark zeromodes profiles, obtained numerically from lattice, to the known analytic solutions within the semi-classical theory of dyons. While these profiles are not in general topologically protected, we still find a surprisingly good agreement between the analytic (semi-classical) and lattice profiles. Continuing forward from our preliminary findings reported in Ref. [28], we describe here our findings in a more extensive form, with several new results.
The main conclusion of our work is that a semi-classical description of the fermionic zero-modes of QCD in terms of an ensemble of dyons is quite accurate, just above T c . In a more broader context, let us note that recent lattice measurements of the topological susceptibility in QCD [29][30][31][32][33] were related to a semi-classical description, in terms of a dilute gas of instantons [34] only at fairly high temperatures, beyond 3 T c . Our study, focusing on a different temperature regime just above T c , sheds light on the instanton ionization phenomenon.
The paper is organized as follows. In section II we describe the basic computational techniques and resources used in this work. Then, in section III we review the known expressions for the corresponding solutions of the Dirac equation (the zero modes) in the dyon background. Next section IV contains our main results, presented in a series of sub-sections. We show our lattice results for the quark zero-modes as a function of temperature for dif-ferent choices of the temporal periodicity phases for the valence fermions. We qualitatively describe the Dirac zero-modes obtained from lattice calculations, emphasizing the evidences that those can be explained in terms of dyons. For overlapping as well as well-separated dyons we show that the lattice fermionic density agrees well with those calculated in the semi-classical theory of dyons. For the first time, we calculate distribution of the dyons in sizes and locations, as a function of temperature. We finally conclude with possible applications of our results for a more comprehensive understanding of the crossover transition of QCD in terms of its topological content.

II. COMPUTATIONAL TECHNIQUES
As already mentioned in the Introduction, our aim is to use lattice fermions with exact chiral symmetry and index theorem on the lattice for understanding the topology in hot QCD. This is achieved in two ways. First of all, the lattice configurations for 2 + 1-flavor QCD at finite temperatures has nearly exact chiral symmetry. These were generated by the HotQCD collaboration using Möbious domain-wall discretization for fermions [36,37] and Iwasaki gauge action, previously used for studies of the chiral crossover transition temperature T c and U A (1) breaking in QCD [3].
The lattices have spatial dimensions N s = 32, with N τ = 8 temporal sites and N 5 = 16 sites along the fifth dimension. For this extent of N 5 , the residual chiral symmetry breaking due to this specific fermion discretization on the lattice, is of the order of ∼ 2 × 10 −3 . The pseudocritical temperature measured using these configurations is T c = 155(9) MeV [3]. We have selected four temperatures T /T c = 1.0, 1.08, 1.1, 1.2 for our study. While these may appear very close to each other, this selection is in fact rather wide, because the relevant quantity for dyons, the average Polyakov loop, varies significantly between near zero and O(1) values. Therefore, our choice is well suited for studies of physical phenomena that depend crucially on the non-trivial holonomy. A more detailed discussion of the holonomy of these analyzed configurations can be found in Sec. IV E.
Some configurations studied were selected to have total topological charge |Q| = 1, and respectively a single zero mode. We found this way is the cleanest set-up to identify and study exactly the zero modes. We have also studied configurations with |Q| = 0, 2, 3 to check the robustness of our results. The details of the configurations used in our analysis is summarized in Table I. Unfortunately, the domain wall fermions used in generating the configurations, turned out to be not precise enough for our purposes. Due to small but finite residual chiral symmetry breaking, the domain wall fermions used in simulation do not satisfy an exact index theorem on the lattice, neither do they have well defined zero-modes. Therefore we have instead used another formulation of  lattice fermions, known as the overlap fermions [38], as the so called valence fermion operator, where D W is the Wilson-Dirac operator for massless quarks in four dimensions and M is the domain wall height which is chosen to be M = 1.8 ∈ [0, 2) to simulate one quark flavor on the lattice. Owing to the fact that the overlap operator has an exact index theorem, even for a finite lattice spacing [39], its zero modes must be related to topological objects resigning in the underlying gauge configurations.
We have implemented the overlap operator on the lattice by explicitly evaluating the sign-function via deflation on the eigen-space of H 2 W consisting of the lowest 25 eigenvectors and approximating on the rest of the eigenspace by a Zolotarev Rational Polynomial with 25 terms. We have ensured that the sign function is realized as precise as ∼ 10 −9 . Furthermore we have checked that for each gauge configuration the resulting overlap-Dirac operator satisfies the Ginsparg-Wilson relation [40], to a precision of 10 −9 or better.

III. FROM CALORONS TO INSTANTON-DYONS AND THEIR DIRAC ZERO MODES
The term caloron is used for the finite temperature classical solutions of gauge fields [41] with unit topological charge and trivial Polyakov loop vacuum expectation value (VEV). A crucial new step made in [11,12] was to generalize it to the case in which the Polyakov loop has a nontrivial VEV. In the relevant case of the SU (3) gauge theory, it can be represented by three phases µ 1,2,3 , called holonomies, such that the average value of the Polyakov loop is P = 1 3 T r[exp[iDiag(µ 1 , µ 2 , µ 3 )]]. These three phases can be measured on lattice configurations, both locally, and as an average over the entire spacetime volume.
The classical gauge field solutions of Yang-Mills equation with arbitrary µ 1,2,3 , based on a combination of Nahm transformation [42] and the Atiyah-Hitchin-Drinfeld-Manin (ADHM) construction [43], is rather involved. Let us just mention that this solution has twelve parameters, out of which nine of them are the spatial locations of the three distinct dyons, with the additional three parameters being the U (1) phases. If the locations of dyons are well separated in space, then for an appropriate gauge choice (the hedgehog gauge), one finds that these dyons are basically monopoles, with the gauge potential A 0 as a scalar field. It is instructive to put these phases µ 1,2,3 on a unit circle, as shown in Figure 1 by small red circles. One of the phases µ 1 of the Polyakov loop is real, and two others are complex conjugates of each other. The explicit values of the phases shown in Figure 1 is a typical representation in the confined phase of QCD. The actions for each of these dyon types are proportional to the lengths of arcs (segments) between two successive phases on the circle, e.g. for M 1 -dyon the action is proportional to µ 2 − µ 1 . Since the three arcs together make a complete circle, the sum of the action of the three individual dyons is always equal to 8π 2 /g 2 , the action of the caloron. Since the solutions are (anti)selfdual, the action and the topological charges are (up to a sign) proportional to each other. Dyons thus have noninteger topological charges.
And yet, from the index theorem, the number of Dirac zero modes can only be an integer. This would naively imply that by counting the number of zero modes of the QCD Dirac operator, we will always get an integer topological charge and can never identify the individual dyons that constitute a caloron. To resolve this dilemma we recall that in QCD at finite temperature, all quark fields acquire a phase φ u = φ d = φ s = π when traversed once along the Euclidean temporal direction. This is the socalled anti-periodic boundary condition due to the Grassman nature of the quarks. However, one can assign any arbitrary periodicity phase φ one wants, to the valence quarks, if we use the valence fermion zero modes merely as a tool, or a filter, to identify the topological content of the gauge configurations. The identification rules are now simple; depending on which segment of the phase circle the φ lies in, the dyon corresponding to that segment produces a normalizable zero-mode of the Dirac operator. For the other dyon backgrounds, the corresponding Dirac operator with phase φ will have nonnormalizable solutions, which are not physical. In this work, we have selected three possible phases for the valence overlap fermions, φ = π, ±π/3 as shown in Figure 1.
Since these values are placed in segments corresponding to L, M 1 , M 2 -type dyons, these three choices would let us identify all three types of dyons via its zero modes. Obviously this method of detecting dyons by tuning the periodicity phase of the valence Dirac operator works best if the dyons are well separated i.e., when the calorons are ionized into its constituent dyons. In this work we will apply this method not only to well-separated dyons, but also to topological clusters, made of considerably overlapping dyons, see Sec. IV.
As it will be clear from the analysis to follow, we will fit the shape of the topological clusters obtained from the lattice with the Dirac zero modes of overlapping dyons derived analytically in the semi-classical approximation. The details of the fitting procedure is explained in Appendix A, where we have numerically calculated the zero-mode density used in the fit, which is outlined in Appendix B. Here we only review the basic theoretical foundations for the derivation of fermion zero modes on dyon backgrounds. It has been derived earlier [44] that the fermion zero-mode density can be related to the Laplacian of a certain Green's function (diagonal part) which is defined on the circle of phases of the Polyakov loop [45]. The fermion zero-mode density ρ(x) at any 4D spacetime point x is derived to be, where the Green's function f x (φ, φ ) itself is the solution of the differential equation, with the corresponding strengths proportional to distances between the centers of the m-th and the (m + 1)-th dyons, which are given as x m , x m+1 respectively where m = 1, 2, 3. The x-dependent potentials are r 2 (x, φ) = r 2 m (x), φ ∈ [µ m , µ m+1 ] defined differently in corresponding segments since r m (x) is the distance between the observation point x and the center of the m-th dyon. Till recently the analytic solution of this equation on a torus was only available for the SU (2) gauge group [44], the case with only two species (L and M ) of dyons. Since here we are working with the SU (3) gauge group, we solve Eq. 4 numerically, details of which are outlined in Appendix B. It is only during finalizing our manuscript we learned that an analytic solution for the gauge potentials in a general SU (N c ) gauge group with a non-trivial holonomy on a torus has been worked out in Ref. [46]. It would be now possible to compute the fermion zero-modes on such backgrounds.
We will henceforth use the numerical solution to fit it to the topological clusters obtained from lattice simulations, and will successfully identify and describe quantitatively different clusters. In particular, the spacetime locations of all dyons, as well as the Polyakov loop values in terms of its phases µ i will be determined in each case. This will enable us to distinguish between overlapping set of all three species of dyons on one hand, and a caloron with a trivial holonomy on the other.
We have implemented the quark periodicity phases φ along the temporal direction for valence overlap Dirac operator by setting the gauge links along the temporal boundary to U 4 (β + 1/T ) → e iφ U 4 (β), corresponding to three different choices of φ = −π/3, π/3, π, the last one corresponding to the usual anti-periodic boundary condition. As evident from Figure 1, these choices of φ ensure that the zero modes are equidistant from each other and thus not too close to each µ m , where its wavefunction density has a discontinuous jump from one dyon sector to the other. For these three choices of φ we have then calculated, at each temperature, the first 6 eigenvalues and eigenvectors of the valence overlap Dirac operator, on the (domain wall) configurations using the Kalkreuter-Simma Ritz algorithm [47]. At the two highest temperatures T /T c = 1.1, 1.2 we have also calculated the first 30 eigenmodes of the QCD Dirac operator to understand the dependence of the eigenvalue spectrum on the quark periodicity phases φ. Two gauge-invariant observables for each eigenstate are its density ρ(x) and the chiral density ρ 5 (x), defined at each spacetime point x as, Studying the spacetime profiles of the zero modes and tracking their change in spatial positions for different boundary conditions on the lattice, we identify the presence of calorons or its substructures, the dyons. The eigenvalue method is numerically unambiguous compared to cooling or Wilson flow techniques, since the lowest eigenmodes are automatically separated from the ultraviolet modes, making them insensitive to large ultraviolet fluctuations in the gauge field. This therefore, allows us to study the infrared modes of the QCD spectrum in a non-invasive manner.
FIG. 2: Density ρ(x, y) of delocalized zero modes for φ = π (red), π/3 (blue) and −π/3 (green) for T /Tc = 1, 1.08, 1.1 at a fixed value of the coordinate z at the maxima of the profile. The Euclidean time coordinate is summed over.

A. Properties of quark zero-modes near Tc
Let us start this section with the following question: are the dyons separate entities, appearing at different locations, or are they simply parts of certain multi-dyon clusters, e.g. the original calorons? In the latter case, the positions of the zero-mode density peaks will remain localized at the same point, whereas in the former case, they will be at different locations, as the valence quark periodicity phase φ is changed.
To answer this question, we look for the zero modes (of the valence overlap operator) for three different temporal periodicity phases φ = π, ± π 3 . We study them on QCD ensembles at 3 different temperatures between T c -1.1 T c . For all the statistically independent configurations considered, we have observed that both types of behavior mentioned above, are in fact present. Exam-ples of such regimes are presented in Figures 2 and 3 respectively. In Figure 2, the peak positions of the quark zero-mode density changes for different temporal periodicity phases φ = π, π 3 , − π 3 providing us the first evidence that dyons are in fact observable as separate independent entities. We observe similar behavior in a significantly large number of configurations at three different temperatures T < 1.2 T c , observing that the probability to find such configurations with distinct shift of zero modes corresponding to well separated dyons, reduces with temperature.
A contrasting picture of the topological cluster is shown in Figure 3. Here peaks of the zero-mode profiles corresponding to different fermion periodicity phases φ do not shift. This may be due to two reasons: i) Firstly, the topological objects are still L and Mdyons, but just located very close to each other. We remind that classically their mutual deformations keeps the total action constant, so any interaction between them should be of quantum nature.
ii) Secondly, it can be that the local value of the Polyakov loop is trivial P ≈ 1 and thus the whole dyon construction collapses back to instantons.
We will unambiguously identify the exact nature of these zero modes in Sec. IV F. The next obvious question to address is the density of different types of dyons, and its dependence on temperature. Comparing ensembles at different temperatures, we observe that at the highest temperature 1.2 T c , the occurrence of fermion zero-modes decreases substantially. Moreover we rarely observe well separated zero mode profiles as we change the quark temporal-periodicity phase φ. This may be due to the fact that at higher temperatures, the holonomy is trivial resulting in the average value of Polyakov loop to be close to unity. One of its eigenvalues ν m = 2π, whereas two others are identically zero due to which only one type of dyon survive at high temperatures in statistically large number of configurations. We further elaborate on this case again in Sec. IV F.

B. Properties of the near-zero modes
So far we have analyzed the most clean case of zero modes in QCD configurations with Q = ±1. Now we proceed to study the near-zero modes in the same configurations. The spectral density of near-zero modes of the Dirac operator is an important observable since it is directly related to the chiral condensate through the Banks-Casher relation. Depletion of their density nearzero eigenvalues thus means disappearance of this condensate and the restoration of the SU A (N f ) chiral symmetry. A separate important topic is the fate of anomalous U A (1) chiral symmetry in this temperature range, also influenced by a certain functional dependence of the near-zero eigenvalue spectrum. The very appearance of Dirac eigenvalues close to zero is due to collectivization of a large number of topological objects, scaling with the lattice volume ∼ O(V 4 ). This phenomenon is related to the density of such objects, but in a subtle indirect way, since topological objects can make neutral clusters with zero overall topological charge. Since near-zero modes arise as a result of interactions of the underlying topological objects, these thus contain crucial information on their nature and the density.
The first qualitative question we address here, is whether the density of different species of dyons are the same or not, as a function of temperature. For this we measure the chiral density profile ρ 5 (x, y) of the first near-zero mode as a function of two spatial coordinates, with the z-coordinate fixed and the temporal coordinate summed over. It is very striking that the densities of L and M -type dyons, observed via peaks in the chiral density profiles of the near-zero modes, are indeed very different! As an example, let us show the near-zero mode ρ 5 -profiles for the configuration at 1.1 T c . We plot the chiral density profiles for two different periodicity phases φ = π, π/3 in the upper and the lower panels of Figure 4, respectively. In the case φ = π we observe only one closely located dyon and an anti-dyon pair, while for the case when φ = π/3, we observe more than one dyon and anti-dyon pairs, both close as well as widely separated.
To explain these results, let us recall that the finite temperature QCD ensembles which are studied here were generated using the usual anti-periodic boundary conditions for the quark fields along the temporal direction. Hence the presence of weakly-interacting widely separated L-dyon pairs are already suppressed due to the fact that in such cases the quark determinant is close to zero. However widely separated M -dyon pairs (which do not correspond to anti-periodic boundary conditions along the temporal direction) can still exist with a small but finite probability, in all such ensembles. When we probe such ensembles with valence overlap fermions with quark periodicity phase π (corresponding to L-dyons) then its near-zero modes arises only due to closely separated pairs of L-dyons and anti-dyons. Changing the valence quark periodicity phase to φ = π/3, its near-zero modes measures the M -dyon pairs which may be either close or widely-separated, all of these cases contribute to the near-zero spectrum. To summarize, the near-zero modes of the Dirac operator are excellent probes of the identity and the interactions between dyons. It allows us to distinguish between L and M -dyon species unambiguously and also to estimate the interactions between them, which we show through this example.

C.
Spectral density of the Dirac eigenvalues with different temporal periodicity phases In the preceding section we have shown that the Ldyons are paired, and experience relatively stronger interaction, compared to the M -dyons. Now we proceed from qualitative examples to a full statistical sampling of all near-zero modes we have, making this comparison more quantitative.
The spectral density of the Dirac eigenvalues, of the valence overlap operator for QCD gauge configuration at 1.1 T c is shown in the top panel of Figure 5. Three sets of points correspond to different quark periodicity phases. Clearly one observes from it that an ensemble of L-dyons, generating near-zero eigenvalues shown in square symbols, is less dense than that by M -dyons (circles and triangles). Furthermore, the corresponding quark condensates (corresponding to the spectral density of, near-zero eigenvalues which extrapolate to zero in the infinite volume limit, ρ λ (λ → 0)) vanishes for L-dyons and remains finite for M -dyons. So, in the real world, in which quarks are fermions and the physical quark condensate is generated due to the L-dyons, this temperature is clearly in the chiral symmetry restored phase. And yet, at the same temperature the ensemble of M -dyons is dense enough to still be below its chiral restoration temperature! The same statement holds for our highest temperature T = 1.2 T c as well, as evident from the eigenvalue spectrum shown in bottom panel of Figure 5.
In order to appreciate how different the ensembles of L and M dyons are, we note that the spectral density is shown up to rather large eigenvalues, λ ∼ 60-80 MeV. In other words, the number of low-lying eigenstates we calculated is ∼ 30 per gauge configuration, a factor 3-5 larger than the mean number of dyons of all kinds in them. This means that the zero-mode zone(ZMZ) generated by such topological solitons is responsible only for a fraction of the spectral density corresponding to roughly ∼ λ/T < 0.2. Moreover, the difference between the spectral densities of L and M -dyons is evidently quite large. Here the eigenvalue density is calculated after performing statistical averaging of more than 30 topologically independent configurations in order to ensure that all allowed topological charge sectors are ergodically sampled. Hence it is a more robust observable insensitive to topological freezing and gauge field fluctuations.
To summarize this section, since at higher temperatures only closely located L dyon-antidyon pairs exist,there is a significant depletion of the corresponding near-zero quark modes. However since the M -dyon pairs of all separations survive, the corresponding fermion near-zero density does not fall with temperature as drastically as that of the L-dyons. We will measure the distribution of the separation lengths between L and M -dyon In this section we provide some evidences that the cutoff effects in our studies are under control and our conclusions in the previous sections are physically robust. To this end, we study the zero and near-zero modes of the valence overlap Dirac operator for different quarkperiodicity phases at ∼ T c and 1.08 T c on the gauge configurations used previously, but now with the highest modes of the gauge fields (momenta ∼ 1/a) removed.
To implement this we performed two steps of HYP smearing [48] on the above mentioned gauge configuration in order to smoothen the ultra-violet fluctuations, which exist on the length scale of one unit of lattice spacing. If the zero modes we observed earlier were due to gauge dislocations i.e., localized artifacts of typical extent of about a lattice spacing, which are known to occur with a finite probability at coarser lattice spacings, then they should disappear after smearing. The topological objects like dyons and calorons are of typical sizes of several lattice spacings (for fine enough lattice spacings they can extend up to several lattice sites) and hence their topological properties should survive such smearing, in fact, any such continuous deformations. Indeed the ultra-violet modes of the gauge fields seems not to affect the shape and location of the zero-modes, as is evident from the comparison of the zero-mode profiles at 1.1 T c shown in Figure 6, before and after minimal smearing and for three different quark temporal-periodicity phases. Furthermore the height, width and the number of zero-modes of the fermions of each periodicity phase φ, remain unchanged after the ultra-violet smearing, which indicates that these are not artifacts due to finite lattice spacing. These conclusions survive even at temperatures close to T c . The near-zero modes too remain unaffected from the smearing of the ultra-violet fluctuations leading us to conclude that the lowest modes of the Dirac operator we use are indeed insensitive to lattice cut-off effects.

E. Multi-parameter fits of the topological clusters
To establish connections between our lattice results for the quark zero-mode density and the corresponding analytic expression in Eq. 3, we calculate the latter numerically by constraining the locations and widths of the dyons from the lattice data. There are several input parameters needed to initiate the fit routine. Among these parameters, the important ones are the holonomies i.e., the Polyakov loop phases µ i , i = 1, 2, 3, since these define the actions (and thus topological charges) as well as physical sizes of the dyons. The details of the fitting procedure are rather technical and are mentioned in Appendix A for the interested readers. Here we present here only the results of direct physical interest.
Among the obvious questions to raise here are whether the values of the holonomies obtained from the fits do show a tendency to vanish at T c (the lowest of our temperature range), and how well the values of the Polyakov loop we get agree with what is known from its continuum estimate in QCD obtained from independent lattice studies. Our results for final values of P obtained as a result of the fits are shown in Figure 7. Agreement between two set of points demonstrate convergence of the fits, starting from very different initial trial values, P = 0 (circles) and P = 1/3 (squares). Furthermore, it displays the expected trend of having a vanishingly small value at T c , growing to P ∼ 1/2 at our highest T = 1.2 T c , although with a large error as the number of configurations (used for the fit) with |Q| = 1 at T = 1.2 T c are quite limited. The final values at each temperature obtained from the fits, agrees quite well with the continuum estimates of the Polyakov loop obtained from independent lattice studies [5,49], shown as a gray band in the same Figure. The other characteristic feature of the topological clusters obtained as a result of the fits, are the distribution FIG. 6: The valence overlap quark zero-mode densities measured on a two-step HYP smeared (red) and the original unsmeared (blue) gauge configuration for three different quark periodicity phases φ = π (top panel), φ = π/3 (middle panel) and φ = −π/3 (lower panel), shows a remarkably good agreement. This is a configuration at 1.1 Tc with a net topological charge Q = 1.
of the dyon locations. With a triad of L, M 1 , M 2 dyons identified, one can measure the distances between the centers of a L and any one of the M -dyons. The histograms of the probability distribution of the distances between these two species of dyons is shown in Figure 8. The most probable distance between a L and M -dyon pair is between 0.2-0.3 fm at 1.1 T c whereas at 1.08 T c it seems to be lower but within errors, are equally probable in the earlier range. Within the statistical errors of our data, we do not observe any significant temperature dependence of the topological cluster sizes. Thus with reasonable consistency the average separation between L and M -dyons in the temperature range studied so far, is about ∼ 0.3 fm.
We also compare the average distances between a L and a M -dyon to that between two M -dyons, results of  which are shown in Figure 9. At all the temperatures above T c , the average distance between M -dyons seems to be always larger but at present the errors on the values are quite large to make very precise predictions. Since the relation between inter-dyon distance r M L and the caloron size parameter is ρ 2 = r M L πT [50], the value, r M L ≈ 0.3 fm obtained from our study is in agreement with mean instanton sizes ρ ∼ 1/3 fm of the Instanton Liquid Model of the QCD vacuum [51].
A few more additional comments regarding the fits are as follows: i) We remind that within a 13-dimensional parameter space of our fit ansatz, it is non-trivial to determine the functional dependence of the χ 2 landscape. It could well be possible for the final set of parameters obtained after the χ 2 minimization, to end up in an entirely different minima. Agreement between the two different fit results implies that, at least on average, we do not get different results for very different initial assumptions of parameters and the fitting procedure is reliable with good convergence property.
ii) We usually did not consider the final fit values from those ensembles where the dyons of the same type are too close to each other e.g. in Figure 9, as the fit results did not yield a low χ 2 .
We conclude the section with a very optimistic observation that within (admittedly large) statistical errors, our fitted values for the average value of the Polyakov loop are quite consistent with the average QCD data at all temperatures studied. Furthermore the fits to the distances between dyons quantify the notion of ionization of instantons, to a typical size of a femtometer. Before we proceed to answer this question, we introduce another observable which distinguishes between a dyon and a caloron profile, allowing us to confirm unequivocally the identity of the zero-modes in our studies as dyons. We specifically calculate the profile of the zeromode density along the x-τ plane, (obtained setting the other spatial coordinates at the maxima of the peak) and propose this as a diagonistic tool to distinguish between dyons and calorons. If we now compare the shape of the density profile obtained from lattice to that calculated from the analytical formula in Eq. 3, both for a caloron (with a trivial holonomy), and a dyon as shown in Figure 10, we find a much better agreement with the latter. In particular, the large distance fall-off of the densities serve as a very discriminating evidence for the presence of a dyon. Near the confining values of the holonomy, fermion density-profile corresponding to a dyon has a faster fall-off while fitting the zero-mode peak accurately, which is consistent with the lattice data. For calorons, however the large distance fall-off is more gradual which clearly shows deviations from the lattice profile. In both these plots the logarithm of the density is shown since it resolves this long-distance fall-off of the profile more accurately. Though the comparison and fit to lattice data is explicitly shown for one specific gauge configuration at 1.08 T c , with a temporal periodicity phase φ = π for the quarks, this trend is generically true for a majority of the configurations we have studied so far, just above T c .
As explained in detail in the Appendices, we have numerically solved for the analytic formulae of fermion zero-mode densities corresponding to mutually noninteracting L and M -dyons as in Eq. 3, matching their peak heights and positions to those of the zero-modes located within the lattice box, for example, the ones showed in Figures 2, 3. The fact that the fits converged with a reasonably good precision, allows us to conclude that a semi-classical theory of weakly-interacting dyons can, quite accurately, explain the QCD vacuum quite well near T c .

V. SUMMARY AND OUTLOOK
This work is inspired by multiple earlier works spread over the last couple of decades, trying to understand the underlying topology of the gauge fields, generated numerically in lattice gauge theory simulations. Determining the exact identity of the topological solitons is a difficult problem, requiring sufficient control on finite volume effects and lattice spacing artifacts. Our work has systematically developed more precise and discriminating tests, allowing us to reveal the identity of the topological objects and clusters present in finite temperature QCD gauge ensembles. We focused on the temperature interval (1 − 1.2)T c just above the chiral crossover transition FIG. 10: The log(ρ(τ )/T 4 ) corresponding to the zero mode of valence overlap fermions, whose profile was shown earlier in the middle panel of Figure 3, for a temporal periodicity phase φ = π. The lattice data for the fermion zero-mode density, shown as triangles is compared to the analytic density profile for a dyon, shown as circles, and for a caloron with trivial holonomy, shown as squares. The filled symbols represent the logarithm of the normalized zero-mode density at a value of the x-coordinate close its maxima whereas the empty symbols correspond to a value of x-coordinate far away from its maximum. The χ 2 310, 423 respectively when the fit to the lattice data is done to a dyon and to a caloron with a trivial holonomy.
temperature, T c . In short, the main conclusion of this work can be summarized by the statement that both individual topological objects, and the topological clusters observed do correspond, with unexpected accuracy, to zero modes associated with the instanton-dyons.
The so-called fermionic method used in this work, allows one to probe the topological lumps and their interactions through the lowest eigenvalues of the QCD Dirac operator, instead of looking directly at (very noisy) gauge fields. This filtering tool has been known for quite sometime [52], and has been used to show that the atoms, i.e. the elementary building blocks of gauge topology, are the dyons [25,26], which have fractional topological charges. While these previous works could detect well-separated dyons, the exact identity of the so-called topological clusters could not be resolved till now.
Methodical improvements we have achieved are based on two pillars. One is the extensive use of overlap fermions which has exact chiral symmetry to very high precision on finite lattices. The other is the use of QCD configurations generated using domain wall fermion discretization with physical quark masses, provided to us by the HotQCD collaboration. It is important that even domain wall discretization was found to be not so accurate for our purposes. To get the purity of chiral and topological properties, we need to apply the overlap fermions, which was a significant change on the algorithm part of the project.
Furthermore, we took a challenging task of devising a fitting procedure, to understand the constituents of topological clusters seen in the lattice simulations. QCD with N c = 3 gauge group has three different species of dyons. With all three dyon species overlapping with each other, the shape of each fermion zero-mode configuration depends on twelve collective coordinates. The lattice data however provided a few thousands of relevant entries, good enough to get stable fits. The fits thus provided very nice explanation of the topological clusters in terms of the overlapping dyons. We found that the corresponding topological structures sometimes consist of three well separated dyons, and sometimes arise from strongly overlapping ones. This qualitative finding is already important, as it rules out a mean-field type picture of the dyon plasma. The dyons are strongly correlated, not randomly placed. The fits turned out to be so good, that we could even reproduce the shape of the zero-modes in terms of a semi-classical theory of dyons. We further checked the robustness of our fitting procedure and found that the average Polyakov loop values constructed out of the final holonomies, extracted from the fit, agrees remarkably well with the corresponding values in QCD. The other advancement we made was a very precise measurement of the near-zero modes. Scrutiny of near-zero modes revealed qualitatively very different interactions between L-M dyons compared to that between M -dyons. Not only do we clearly see that there are more M -type dyons than L-type ones, our results for the Dirac spectral density ( Figure 5) show that while the latter form finite neutral clusters, i.e., closely spaced LL pairs and have zero quark condensate, the M -dyons still do break chiral symmetry up to 1.2 T c .
The conclusion that emerges out of our study is highly nontrivial: the topological lumps we study here, are located within a quark-gluon plasma, at temperatures above T c , with a large density of thermal quarks and gluons around. There is no doubt that the gauge fields of these lumps are not close to classical solutions (minima of their action) but strongly deformed by plasma excitations. And yet, the corresponding Dirac zero eigenstates, in number and even exact shape, appear completely unperturbed by them, in accurate correspondence to what one gets from solving the Dirac equation for weakly-overlapping classical dyons. This conclusion is also practically relevant; if this phenomenon survives further scrutiny near the continuum limit, it would imply that one can use semi-classical ensembles without thermal corrections, at least in the fermionic sector, to study topological properties in QCD.
Needless to say, one can improve the statistics and extend the temperature range of this study, finding in particular, the dyon densities at all temperatures, and elucidating their exact role in the formation of quark condensates. One can also deform QCD by using nontrivial flavor-dependent periodicity phases φ f (or using imaginary chemical potentials) in simulations and study the intimate connection between deconfinement and chiral phase transitions. Another challenge is to get a better understanding of the effective interactions between the dyons. Since dyons have magnetic charges, one may wonder if the magnetic confinement, known to be present above T c , have a role in it. In other words, it will be interesting to investigate how dyon-ensembles at high temperatures may interpolate to Polyakov's monopole-driven confinement in three dimensional gauge theories.

VI. ACKNOWLEDGMENTS
This work was supported in part by the Office of Science, U.S. Department of Energy, under Contract Nos. DE-FG-88ER40388 (E.S) and de-sc0012704 (R.N.L) and by the Department of Science and Technology, Govt. of India through a Ramanujan fellowship (S.S). We thank the HotQCD collaboration, formerly also consisting of members from the RBC-LLNL collaboration, for sharing the domain wall configurations with us. S.S. would like to thank Margarita Garcia Perez and Antonio Gonzalez-Arroyo for many insightful comments and helpful suggestions on the first version of the draft. We gratefully acknowledge the computational facilities at the Institute of Mathematical Sciences. Our GPU code is in part based on some of the publicly available QUDA libraries [53].
Appendix A: Details of the fitting procedure In this section, we explain the details of how we have numerically fitted the solution of a dyon density profile obtained from Eq. 3 to the density of the fermion zero eigenmodes, measured on gauge configurations on the lattice. The details of how we obtain the analytic solutions corresponding to Eq. 3 are given in the Appendix B. We chose the (overlap) fermion zero-mode solutions on gauge configurations with a net topological charge |Q| = 1, for the fit. The fit was performed simply through a gradient flow minimization of the χ 2 function, made possible thanks to a large number of data points available from our lattice calculations. However in this procedure, there is a finite probability that the fit procedure may result in ending up in a local minimum of the χ 2 landscape as a function of the fit parameters. In order to circumvent this possibility, we made an initial guess for a couple of parameters based on the properties of the zero mode. For example, we set the temporal periodicity phase and the spatial co-ordinates of the dyon to be closest to the maximum or the peak of the fermion zero-mode density under study, such that the initial shape of the analytic solution mimics quite well the overlap zero-mode profile. In Figure 11, we show typical zero-mode density profiles corresponding to a L-dyon as a function of any two space coordinates, which is simply a numerical solution of Eq. 3 with the two other M-dyons either placed close or far from each other as well as the L-dyon. As evident from the Figure,  profile as close to the realistic profiles obtained from lattice already improves the final fit.
In total we had thirteen parameters for the fit, nine of them being the spatial coordinates of three dyons respectively, two of them the free angles or the holonomies µ i , the temporal periodicity phase and the normalization (height) parameter. Even the spacetime profile for the zero mode on |Q| = 1 configurations does not always comprise of a single distinct peak but may have additional smaller peaks often far enough away from each other. Thus for the analytic fermion zero-mode solution to be a good fit, we therefore considered the height of the peak as a fit parameter that could take any value between zero and one, to take into account that the peak only caries a fraction of the total normalization. The error in estimating the zero-mode density at each spacetime point on the lattice was assigned to be identically equal, which implies that the points near the maxima of the density profiles are relatively important in deciding the convergence of the fits, since all solutions approach zero values at large distances, albeit with a different dependence on the coordinates. We found that the initial choice, in which we construct all three dyons starting from Eq. (3), close to each other such that the relative separation is ∼ 0.05 fm (like shown in Figure 11), is the most optimal. This is due to the fact that during χ 2 minimization, it is more likely for the dyons to move further away from each other. Furthermore, the fits were performed with two different initial choices of average Polyakov loop P = 0, which correspond to the choice of the holonomies µ 1 = 0, µ 2 = 2π/3, µ 3 = 4π/3 and P = 1/3 for which the holonomies are µ 1 = 0, µ 2 = π/2, µ 3 = 3π/2. The choice of P = 0 is closest to the corresponding values physical values just above T c , however P = 1/3 was chosen as an additional initial guess to check how our fit procedure performed when we are sufficiently far way from the anticipated physical value. It was however ensured that the values µ i did not cross from π/3 and −π/3 corresponding to two distinct dyon solutions. As evident from Figure 7 the final results for the average Polyakov loop agree within errors, for either choice of the initial values, giving us a confidence on the robustness of our fitting procedure.
We have also checked the sensitivity of our numerical fitting procedure on different initial choices of the positions of dyons. The final average values of the Polyakov loop, obtained as a result of the fit, showed little dependence on the initial separation of dyons. However the final average distance between dyons resulting from the convergence of the fit did show a mild dependence on the initial separations. For different choices of the initial positions, we obtained the final separations in the range 0.3-0.4 fm. This is related to the fact that a fermion zeromode solution for a particular dyon is most sensitive to position of the dyon situated furthest away from it. It was however assuring that the final χ 2 values were not very different, giving us a confidence in our numerical procedure. This may be as well because the third dyon which is situated anywhere in between, does not affect the outcome of the final fit. While we did derive an analytic solution to Eq. 4 for the dyon-density, the final expression turns out to be quite complicated due to the fact that eight independent boundary conditions need to be fulfilled. Hence substituting the values of the initial holonomies µ i and the coordinates in that expression and using it as the initial trial solution of our numerical fitting algorithm was much more cumbersome. Instead, it is computationally much faster to solve the density from Eq. 4 using standard numerical techniques. For this, we first solve the homogeneous part of the equation (away from the location of the delta functions), which gives us solutions of the form, f x (φ) = C 1 e −(r−iτ )φ + C 2 e (r+iτ )φ where C 1 and C 2 are spacetime dependent functions but independent of φ. There are four such solutions corresponding to the four regions in between the delta functions, hence we have to calculate eight independent parameters C 1 , · · · , C 8 , which satisfy, where c is a column vector that contains C 1 to C 8 as its entries. The operator on the L.H.S of Eq. 4 acting on f x can be represented as a 8 × 8 matrix M (x), containing information about the positions of the dyon zero-modes, acting on the column-vector c. This matrix-equation can be easily solved using linear algebra techniques giving us the solution for vector c. One just needs to calculate the inverse of matrix M which is computationally inexpensive since it is only of size 8×8. The spacetime derivatives of the vector c hence follows immediately as, We recall that, in order to calculate the zero-mode density from Eq. 3, it is sufficient to calculate up to second derivative of the functions C i with respect to the spacetime coordinate x. The computational cost for calculating the density is mainly due to the matrix inversions, since matrix multiplications with the derivatives of M are very fast. The speedup in this process of computing the zero-mode density is around 36 times as compared to using the analytic solution to Eqs 3,4. This method is also highly parallizable, since the dyon zero-mode density at one position is independent of the density at any other spacetime point. We ran this code on four threads of a single processor using openMP, but it can be easily scaled to many more CPU cores.