QCD Topology at Finite Temperature:
Statistical Mechanics of Selfdual Dyons
Abstract
Topological phenomena in gauge theories have long been recognized as the driving force for chiral symmetry breaking and confinement. These phenomena can be conveniently investigated in the semiclassical picture, in which the topological charge is entirely carried by (anti)selfdual gauge configurations. In such an approach, it has been shown that near the critical temperature, the nonzero expectation value of the Polyakov loop (holonomy) triggers the “Higgsing” of the color group, generating the splitting of instantons into selfdual dyons. A number of lattice simulations have provided some evidence for such dyons, and traced their relation with specific observables, such as the Dirac eigenvalue spectrum. In this work, we formulate a model, based on oneloop partition function and including Coulomb interaction, screening and fermion zee modes. We then perform the first numerical Monte Carlo simulations of a statistical ensemble of selfdual dyons, as a function of their density, quark mass and the number of flavors. We study different dyonic twopoint correlation functions and we compute the Dirac spectrum, as a function of the ensemble diluteness and the number of quark flavors.
I Introduction
Topological phenomena in gauge theories have been discovered more than three decades ago, and remain the subject of intense theoretical research ever since. In particular, magnetic objects (monopoles) have been identified as a possible source of confinement Mandelstam:1974pi ; 'tHooft:1981ht , while instantons have been proposed as the driving mechanism for chiral symmetry breaking Shuryak:1981ff ; Diakonov:1985eg .
The index theorem establishes a direct connection between the vacuum topology, and zeroeigenvalue solutions of the Dirac equation, i.e. the socalled fermionic zeromodes. These quark states are insensitive to any perturbative fluctuation of the gauge field, hence encode purely nonperturbative QCD dynamics. Furthermore, lattice simulations have shown that the Dirac eigenstates with nearzero eigenvalues — also known as the ”zeromode zone” (ZMZ)— directly correlate with local fluctuations of the topological charge density. After filtering out quantum fluctuations, lattice fields reveal nearly (anti) selfdual smooth fields responsible for topology and ZMZ states Gattringer . Using only fermionic states attributed to the ZMZ (a tiny subset of Dirac eigenstates, of only about of all eigenstates) one finds the correct pion mass, quark condensate as well as many other hadronic properties. On the other hand, filtering out the ZMZ states removes the chiral symmetry breaking and leads to drastic changes in the hadronic spectrum computed on the lattice. In particular, some masses get shifted by as much as and paritydoublets appear (for a recent analysis, see e.g. Glozman:2012hw ).
This body of results coherently support a picture in which the nonperturbative chiral dynamics in vacuum is mediated by instantons. Indeed, instanton model calculations (for a review see Schafer:1996wv ), have been very successful in reproducing the mass and electromagnetic structure of pions pionFF , vector mesons hadronicresonances , nucleons Negele ; diquarks ; masses ; nucleonFF1 ; nucleonFF2 ; nucleonFF3 and even the rule for hyperon delta12 and kaon delta12kaon nonleptonic decays.
In the instanton picture, the width of the ZMZ depends on the size of the typical “hopping” matrix element of the Dirac operator between two instantons, which is of the order 20 MeV, where is the typical instanton size and the typical interinstanton density Shuryak:1981ff . This value is comparable to the typical light quark masses used in may lattice simulations, and this explains why the corresponding results display significant deviations from the naive chiral perturbation theory predictions. Furthermore, the specific shape of the density of eigenvalues in the ZMZ depends crucially on the theory parameters, such as the number of light fermions .
In this work, we will further investigate topological phenomena in the semiclassical picture, focusing on temperatures close to those at which the expectation value of the Polyakov line
(1) 
drastically changes from 1 to 0.
The gauge invariant expectation value (1) defines the holonomy of the gauge connection corresponding to a full circle around the periodic time direction and is related to the freeenergy of a single static quark:
(2) 
Hence, a drastic suppression of reflects the onset of the confinement phenomenon and we shall denote the corresponding critical temperature with .
From Eq. (1) it follows that a nontrivial holonomy (i.e. ) reflects a nonvanishing vacuum expectation value (VEV) of the component of the gauge field. In the semiclassical picture, this condition is not fulfilled by the standard instantons and calorons (finite temperature instantons). Indeed, in such topological classical solutions, the component of the gauge field vanishes at spatial infinity, .
Classical gauge configurations, generalizing the instanton solution to a nontrivial holonomy case, are called the KraanvanBaalLeeLu (KvBLL) calorons Lee:1997vp ; Kraan:1998sn . These authors have shown that instantons basically split into certain substructures, the “constituent dyons”.
The collective variables for the constituent dyons are different from that of standard instantons: while each instanton is identified its position, size and color orientation, a dyon is specified only by its position and its Abelian phases. The solution and the measure for the KvBLL configurations correctly returns to the standard instanton one in the limit of vanishing holonomy. For review of instantondyons see Bruckmann:2003yq , Diakonov:2009ln and references therein.
In general case, the holonomy can be parametrized in terms the VEV of at spatial infinity, according to the usual Cartan subalgebra notation. Such an asymptotic field can be taken to be spatially constant and diagonal in color space,
(3) 
The latter condition follows from imposing a null trace.
Inserting such a constant field plus and its quantum fluctuations into the gauge action, one finds that the YangMills commutator generates mass terms for all nondiagonal components of the gluon field. On the other hand, diagonal components of the quantum gauge field commute with the holonomy and remain massless.
It is convenient to adopt a cyclicsymmetric notations such that and introduce the differences between subsequent holonomy parameters,
(4) 
These differences determine the masses of the nondiagonal components of the gauge field. Since the dyon cores are made of such field components, the ith dyon has the core size . On the other hand, the longrange dyon interaction is due to massless diagonal fields and is governed by their electric and magnetic charges (see e.g. Ref. Diakonov:2009ln ).
In this paper, we restrict our attention to the SU(2) gauge group. In this case, the holonomy parameterization simplifies, since there is only one diagonal generator and a single holonomy parameter, . In Fig.1 we show the locations of the holonomy eigenvalues at two different temperatures. At high (case (a)) the holonomy is close to a zero phase, so is small. Respectively, the mass of the first dyon, called M, is also small. The mass of the second dyon is proportional to the complementary part of the phase circle , which at high takes nearly the whole circle. The corresponding dyon is thus heavy.
As , the holonomy moves toward of the circle, as shown in case (b): and in this case both dyons have the same mass. This position corresponds to phases and the Polyakov line vanishes , which physically corresponds to the confinement of heavy charges (quarks).
In Table I we list the quantum numbers of the possible dyon solutions for the SU(2) gauge group. In such a gauge theory there are 4 distinct dyons, denoted with , , and , which are characterized by different (color) magnetic and electric charges. For the general case there are socalled “static” dyons with all diagonal charges and just one socalled “twisted” dyon .
From the technical point of view, the statistical physics of the dyon ensemble is quite challenging compared to that of the “instanton liquid” Schafer:1996wv . This is mostly because, on top of topological phenomena induced by the zeromodes of light fermions, one has to deal with longrange Coulomb forces and even linearly growing potentials due to a “screening” effect, as will be discussed below.
A first, qualitative discussion of dyonic ensemble has been performed by T. Sulejmanpasic and one of us Shuryak:2012aa . It was pointed out that a dyonic system can realize several distinct phases, ranging from a gas of individual clusters at high , to a stronglycoupled liquid plasma at . Ref. Shuryak:2012aa also includes an extensive discussion of the latticebased phenomenology related to selfdual dyons. Extension of that work to the case of adjoint fermions will be given in Tin2 .
The present paper can be considered as a quantitative continuation of that study, we shall define a partition function of the dyonic ensemble which is amenable to practical manybody computer simulations, perform the first numerical simulations and report on the results.
name  E  M  mass 

+  +  
+    
    
  + 
The partition function of the dyon ensemble can be schematically written as
(5) 
where the first factor is the product of the differentials of all the collective variables, the exponent contains the classical actions of all dyons, while the two subsequent determinants of matrices and are related to the bosonic and fermionic zeromodes, respectively. The two remaining (primed) determinants correspond to the nonzeromodes of the oneloop bosonic and fermionic operators, describing the small linearized perturbations of the classical solution. (As usual, the latter determinants are divergent and need to be regularized, e.g. in the PauliVillars scheme, so that the UV bare coupling in the action is combined with the UV cutoff into a physical renormalized charge.) A quantitative characterization of the dyonic ensemble can be achieve by performing Monte Carlo sampling of the partition function (5).
Before we come to such simulations, let us outline the physical phenomena to be addressed. As the temperature drops from to , the value of the holonomy changes from zero to its “confining value” . In the former case the action of the dyons is small and of the ones is large, while at their actions become nearly equal. At high “heavy” dyons can be bound with antidyons by the fermion exchanges into certain clusters Shuryak:2012aa , with zero magnetic but double electric Abelian charge. The formation of such clusters are reflected in the doublebump shape of the Dirac eigenvalue spectrum and a gap around zero: in the gas of clusters the chiral symmetry remain unbroken. As decreases, the action per dyon decreases and their density grow. Also the mass of the fermions decreases. Eventually, at some , the density is high enough to bridge the gap in the Dirac eigenvalue spectrum and produce finite density of eigenvalues at zero. This means that the chiral symmetry gets spontaneously broken.
Let us now briefly outline the structure of the paper. The next section describes bosonic interactions of the dyons, which are separated into approximations made to describe the so called moduli space metric associated with bosonic zeromodes, and electric screening associated with the nonzero bosonic modes. The next section describes our approximation for the third essential ingredient of the partition function, the fermionic determinant.
The subsequent two sections provide further technical details about the partition function (5). Some introduction about the instanton dyons and the setting used in our simulations is reported in section IV. The results of the simulations are presented and discussed in section V, while the main conclusions are summarized in VI.
Ii Bosonic dyon interactions
ii.1 Moduli space metric
In the selfdual (or antiselfdual) subsectors, the dyons are socalled BogomolnyPrasadSommerfeld (BPS) protected objects, which means that their total action is completely determined by the total topological charge. Thus, at the level of classical solutions of the YangMills equations, the total action is independent of the dyon collective coordinates. The corresponding degenerate manyfold of all classical solutions with a particular topological charge , known as the “moduli space”, has certain metric, which can be calculated from the bosonic zeromodes of the classical solution in a standard way. The resulting geometry is not simple: for monopoles it was derived in classic works of Atiyah and Hitchin AH . Even a configurations consisting of two monopoles/dyons gives raise to a quite nontrivial manifold, with a specific singularity known as “bolt” (or arrowhead). Manybody moduli are extremely complicated, and to our knowledge their exact metric and geometry has not been explored in any detail.
Obviously, one needs to make some approximations, in order to perform practical simulations. Following Gibbons and Manton GM and Diakonov Diakonov:2009ln , the invariant volume element for moduli metric can be approximated by the determinant of a certain matrix which is defined in order to incorporate the known limiting cases.
(6) 
The Jacobian determinant of a single LM dyon pair has been calculated by Diakonov, Gromov, Petrov and Slizovskiy (DGPS) Diakonov:2004jn . The nonzeromodes lead to the so called screening phenomenon, which we discuss in section II.2. For the SU(2) gauge group, the matrix is
(7) 
where the matrix raws and columns correspond to the two dyon types and . We emphasize that the expression (7) turns out to be correct at all distances, not only for asymptotically large separations. It was also shown that in the limit of trivial holonomy (i.e. for ) or at vanishing temperature, the resulting measure with (7) reduces to the standard ’t Hooft singleinstanton measure.
For sake of completeness we report also the generalization of Eq. (7) to an arbitrary group,
(8)  
where the indexes run over the different types of dyons.
The effective potential generated by the moduli metric can be defined as usual as . Upon expanding it in powers of , one recovers the usual Coulombic potential at large distances.
A classical Coulomb gas is known to be unstable: particles can fall onto each other. Fortunately, a pair of nonidentical dyons is not affected by such a problem. Indeed, the Coulomb term is complemented by the higher powers of of opposite signs. As a result, the total effective potential turns out to be only logarithmically divergent at small distances, thus the integral which appears in the partition function,
(9) 
is well convergent. This expression implies that the probability of finding two dyons at the same point vanishes.
Let us now discuss the case in which there are dyons of the same kind. Simplifying the famous AtiyahHitchin manymonopole metrics AH , Gibbons and Manton GM have shown that the weight for identical dyons reads:
(10) 
where
(11) 
Note the inversion of the signs here, as compared to the case of nonidentical dyons. We emphasize the this expression is exact only for large interdyon separations, and corrections may be important at small distances.
For a pair of identical dyons, , which remains positive only for distances larger than some “core size” . Note also that the volume element vanishes at the core, which follows from the famous “bolt” geometry.
Diakonov Diakonov:2009ln nicely combined the metric tensors for differentkind and samekind dyons into one symmetriclooking metric, describing the moduli space of arbitrary number of dyons. Let their numbers are dyons of kind 1, dyons of kind 2, dyons of kind . The result is a matrix whose dimension is the total number of dyons and reads:
(12) 
In this expression denotes the coordinate of the i dyon, of kind .
Let us now discuss the shortdistance behavior. Now, the effective potential has multiparticle terms with all powers of . To make sense of it, it is instructive to first calculate it for some examples of configurations. For example, we can consider a “square” consisting of 2 and 2 dyons. In is found that, in all cases, their combined effect can be described by a weakening of the Coulomb divergence.
When the distances between dyons become too small, the Coulombic term becomes unreasonably large, and it needs to be regulated: after all, a dyon is not a point charge and has a certain size. In the present work, we adopt a regulation based on the substitution
(13) 
where we have chosen in the simulation the value of the cutoff parameter , leading to vanishing measure for two identical dyons, in the limit (to be referred to as “soft core” ). This is in contrast with imposing that the measure should vanish at some finite , as in the famous AttyahHitchin metric (to be referred to as “hard core”).
The configuration space of a multidyon ensemble can be sampled with a standard Metropolisbased Monte Carlo algorithm. Within such an approach, we shall not only enforce the global positivity of the metric determinant, but also we will require each eigenvalue of to be positive. This condition is becoming restrictive at high dyon densities, i.e. when fraction of the total volume filled by dyons is , eventually making the Monte Carlo sampling of a very dense dyon ensemble rather inefficient. (This kind of computational difficulty is well known, e.g. in simulations of classical liquids with hard repulsive cores.)
ii.2 Nonzeromodes and electric screening
The DGPS determinant of the nonzeromodes Diakonov:2004jn for a single LM dyon pair contains the socalled “confining” term, proportional to the interdyon separation and to the “screening mass”
(14) 
which displays the socalled electric Debye mass , in the SU(2) gluondynamics. Upon generalizing to any and number of fundamental fermions this expression becomes Shuryak:1977ut .
(15) 
We now briefly elucidate the origin of the linear potential which appears in Eq. (14). Let us start with the quartic term of the YangMills Lagrangian, and imagine two of those fields are classical ones, say from some dyon solution, and two others belong to thermal perturbative gluons of the heatbath. While the latter produce the Debye mass, the former generate the linear dependence on , from the volume average of the squared Coulomb largedistance potentials
(16) 
Note that we consider here neutral pair (i.e. the splitting of one instanton), as a result of which their Coulomb potential cancels out at large distances from the pair and the volume integral is convergent. Clearly, nonneutral configurations cannot be treated this way.
We further note that the form (14) can be obtained directly by the instanton screening term calculated by Pisarski and Yaffe Pisarski:1980md by recalling that the instanton size and the separation are related by the expression
(17) 
which relates the “4d dipole” of the instanton field to the “3d dipole” of the dyon pair made of opposite charges.
Let us now work out the corresponding general formula for screening potential which holds in the manybody case. The sum over all dyonic contributions to can be written as
(18) 
where now the sum runs over all dyons with is the charge and .
One can write as a double sum, in which we separate the diagonal terms from all nondiagonal terms. Unless total neutrality is ensured
(19) 
this integral is divergent at large distances. If there is an overall neutrality, one can regulate the sum term by term, by subtracting the corresponding (rindependent) divergency. Let us say discuss a nondiagonal term which we can then be rewritten as
(20) 
reducing it to the 2body case above (16). Therefore the total answer is a simple generalization of the linear potential
(21) 
with being the electric charges. In contrast to linear potential induced by the confining flux tubes, in this case there is a sum over all pairs of dyons. The sum therefore has terms for dyons, of various signs, which together enforce local neutrality of the dyon ensemble, making as small as possible. It does not diverge in infinite volume because distant dyons may be combined with their counterparts of the opposite charge, producing dipole fields. On a manifold with the topology of a sphere extra care is needed, as discussed in Appendix B.
Iii The fermionic determinant
In evaluating the fermionic determinant, the Dirac operator is approximated by retaining only the contribution evaluated on the subspace of fermionic zeromodes of the individual pseudoparticles :
(22) 
where
(23) 
This scheme was well tested in the framework of the instanton liquid model, where it corresponds to summing up all loop diagrams created by ’t Hooft effective Lagrangian. One of the most important lessons from that was existence of rather narrow “zeromode zone” (ZMZ) of the Dirac eigenvalues 20 MeV, which drive all the effects related to chiral symmetry breaking. The narrowness followed from relative diluteness of the instanton ensemble.
The dyonic formulation also generates a narrow ZMZ: in fact, at nonzero holonomy, all light quarks become effectively massive, so the so called “hopping amplitudes” – i.e. the nondiagonal matrix elements of the Dirac operator between different dyons – are exponentially suppressed with distance,
(24) 
where is the interdyon distance. is the socalled holonomy mass , should not be confused with the ordinary (i.e. current) quark mass . In particular, the former does not break chiral symmetry, since it arises from the interaction of the quarks with the component of the gluon field. An estimate of the magnitude of for various fermions followed from the zeromode solution can be found e.g. in Ref Shuryak:2012aa . Antiperiodic fermions in SU(2) have zero modes on dyons, with the (holonomyinduced chirally symmetric) mass given by
(25) 
proportional to the distance from the fermion boundary condition angle ( the (yellow) lighter circles on the left at the phase in Fig.1) to the (blue) dark circles indicated the holonomy values. If the fermionic boundary phase is changed, the yellow circle moves and then one should select the smallest distance between the two holonomy values.
The most important qualitative trend is the following: as the temperature decreases from high to , the holonomy changes as indicated in Fig.1. As a result, the effective fermion mass decreases and the “hopping amplitude” (24) increases. Another qualitative tendency is of course the reduction of the effective coupling in the dyonic action, which additionaly increases the dyon density and reduces the interdyon distance . Both effects eventually conspire towards breaking the chiral symmetry.
Some of the (block) entries in Eq. (23) are identically vanishing. In particular, the antiperiodic conditions on fermions imply that zeromode exist only for and type dyons, see e.g. Shuryak:2012aa for explicit solution. Hence, we can ignore the overlap zones involving and type dyons. In addition, the Dirac operator mixes only modes with opposite chirality, hence also . Hence, the Dirac operator reduces to:
(26) 
where the the element in block is given by the (approximate) formula Shuryak:2012aa containing the “holonomy mass”
(27) 
Notice that the normalization constant is irrelevant for Monte Carlo ensemble simulations: it may only be needed if some physical significance will be ascribed to the quark mass rather than pushing it to zero. For massless fermion flavors the determinant simply reads
(28) 
iii.1 The case of a single dyon molecule
Let us first consider the case in which there is only one dyonic ”molecule” made by , ,  and  type dyons. The fermionic determinant can be viewed as the secondorder oneloop diagram in ’t Hooft effective Lagrangian, with two vertices with fermionic propagators in between. From Eq. (28) and (27) we get
(29) 
Neglecting the logarithmic dependence, the formula for the determinant for 1 molecule reads as en effective potential
(30) 
where again we have dropped an irrelevant normalization factor. Thus, fermion exchange also creates a linear confining potential, but now between dyons and antidyons. Also, note that there is no minus sign and that an additional parameter – the number of flavors – appears as a factor. Therefore, one expects that fermion exchange to generate the tightly bound clusters first described in RefShuryak:2012aa .
iii.2 The case of two dyon molecules
This is the simplest example demonstrating how dyonic clusters undergo mutual repulsion. To see this we study the two dyon molecule case, restricting to the Ldyon zone and labeling the pseudo particles with , , and we find:
(31) 
Hence
(32) 
from which one immediately sees that the configuration weight vanishes when the two molecules overlap. In the sext sections we shall we shall perform numerical simulations which will allow us to investigate the consequences of such a “Fermi repulsion”.
Iv The setting
iv.1 The setting and parameters
Let us start by enumerating parameters of our physics problem, which we will group into global, internal and external ones.

The global parameters are the number of colors of the theory and the number of fundamental fermions of degenerate current mass . In this work and . One may in addition consider other types of fermions, e.g. adjoint ones, but we leave such an extension to future work.

The internal parameters are those on which the partition function depends, such as the heatbath temperature , and the flavor chemical potentials . Although the theory allows for finite density Monte Carlo simulations, for now we will only discuss the zero density case, , using as the single internal parameter.

The external parameters are those which may be imposed. For example, one may consider modifications of the QCD vacuum in the nonzero (QED) magnetic field: it is now in fact a quite active research field. Another important example of an external parameter is the CPodd angle , which however have a sign problem and need reweighing or analytic continuation. None of it will be discussed in this work.
The dyonic model is defined by the partition function defined above and can be studied for any values of such parameters. After this is done, the results can be mapped/compared to lattice data using appropriate values of these parameters as inferred from lattice simulations at the corresponding temperature.
One such input is the value of the holonomy parameter at temperature . There are rather accurate lattice data defining it, see Appendix C.
Other simulation inputs are the densities corresponding to dyons of various kind . We will also use their proper volumes defined by their inverse densities, i.e. . At the time of this writing the lattice studies of the dyons are too few to allow for a systematic mapping, of the dyon densities at different temperature. While we do not really use such data, for the orientation of the reader we put brief review of the information we find in the Appendix D.
We will only mention here that since the masses of and dyons are different, their densities on the lattice are different as well. Note that such disparity between the densities of the and dyons produces a nonzero density of the Abelian electric (but not magnetic) charge. This is of course not a problem in general, as it is compensated by an asymmetry of the densities of the electrically charged gluons.
It would however create a technical problem for our simulation setting, as we focus entirely on the topological (dyonic) sector. Furthermore, one has to use a compact manyfold – a torus or a sphere – in which all the charges must be compensated and field lines accounted for: hence only totally neutral systems are allowed. One possible solution to this problem is simply adding to each charge the homogeneously distributed compensating charge: see Appendix B for details. Another practical solution we adopt in this work is to ignore the extra dyons: so we simulate an equal number of the dyons.
Finally, about our units. Standard nonperturbative normalization of the gauge field eliminate the coupling from the YangMills equations. It appears in classical action , but since its total value depends on the total number of dyons which is not changed in our configurations, it is absorbed in the overall input density. The oneloop bosonic and fermionic determinants we actually simulate do not contain the coupling. Thus, there is no running and scale dependence in the model at all, and all quantities are essentially dimensionless. Or, as we do, all quantities can always be made dimensionless by multiplying by appropriate powers of temperature .
In particular, the simulations are done for fixed number of the dyons on a sphere, with its radius characterized by the (dimensionless) radius , or by the proper volume per dyon in Tunits, , the inverse of the diluteness (we remind that the total volume is . The values for which the simulations has been done are listed in the Table 2. Note that the density or volume per dyon changes by about two orders of magnitude. The largest values of the proper volume (per dyon) correspond in physical ensemble to weaker coupling, larger dyon masses and higher temperatures. Going down in this volume one goes into the denser systems, stronger coupling and the lower . Respectively we will monitor in simulations how all observables depend on the proper volume, focusing on the correlations between the dyons and the dirac eigenvalues spectrum.
dyons  

64  4.5  28. 
64  3.0  8.3 
64  2.5  4.8 
64  2.2  3.28 
64  1.5  1.04 
64  1.2  0.53 
64  1.  0.31 
The only dimensional parameter which enters the fermionic determinant is the current quark mass , which in our units corresponds to the ratio . Keeping its constant means a temperaturedependent mass. This does not represent a problem as long as one is interested in the massless limit, , in which the classical action becomes scale invariant. Yet in reality, following a prescription which is commonly adopted by lattice practitioners, the actual value of this mass used in our simulations is not zero and is not directly related to physical light quark masses: its role is to prevent influence of the finite volume effects. In particular, we shall will see below that, in order to correctly monitor the chiral symmetry breaking, we are bound to consider masses
V Simulation results
One standard observable in the simulation is the acceptance probability of Metropolis steps, which should be tuned to being comparable to rejection probability to run it efficiently. We estimate the number of Metropolis steps needed to achieve thermalization by standard autocorrelation analysis based on the dyons’ position. Our typical runs include a total of 256 independent configurations, obtained from 32 independent Markov chains, each consisting of 64000 Metropolis steps. All simulations were performed at the Interdisciplinary Laboratory for Computational Science (LISC), at Trento.
v.1 The spatial correlations between dyons
One of the benefits of going from 4d torus to the geometry is that in the latter case it is much easier to define the interparticle distance
(33) 
where is the angle between two points defined via the scalar product of their unit position vectors on and is its radius.
We collected histograms of the various twoparticle correlation functions as a function of the (dimensionless) distance on the sphere are normalized by dividing out the volume element corresponding to random occupation on the sphere. A sample of the results for different values of and different dyonic densities is shown in Fig. 2. One can immediately identify a number of feature which are common to all such correlation functions:

repulsive correlation between identical dyons,

attractive interaction between the dyon and antidyon;

the radius of the correlation decreases as grows
v.2 The Dirac eigenvalues and the final size effects
Our present computer resources allowed us to perform simulations with 64 dyons. In order to estimate the magnitude of finitevolume effects, we performed also a few quenched runs with ensembles consisting of 128 dyons. Two distributions of the Dirac eigenvalues for 64 and 128 dyons performed at the same density conditions are compared in Fig.3. One can see that the finite volume effects are significant only for the smallest eigenvalues. More specifically, one finds distortions of the spectrum for the first bins, say for the shaded area
(34) 
Notice also that that larger systems have smaller finite size effects, as expected.
A rather crude account of the finitevolume effects can be thus made by simply disregarding the data points in the strip identified by the inequality (34), yet we will need more refined fits when we introduce a nonzero quark mass in the next section.
v.3 Quark condensate dependence on the current quark mass
Let us start reminding that there are two types of quark masses: “current mass” generated by the interaction with the Standard Model scalar (Higgs) VEV with the leftright structure
(35) 
and thus violates the chiral symmetry explicitly. The other “mass” is generated by the holonomy. As it is due to the VEV of the vector field component , its structure is chirally diagonal (symmetric)
(36) 
We recall that appears in the parametrization of matrix elements of the matrix, which defines the Dirac operator in the background holonomy potential. (For clarity, we do not use the notation of “constituent quark mass” induced by chiral symmetry breaking in this paper.)
The current quark mass is important for our study because it “seeds” the chiral symmetry breaking. As it is well known to lattice practitioners, in order to study this spontaneously broken symmetry one formally should take the double limit, the infinite volume and massless quarks , in this order. In practice, all simulations are performed on a finite volume, and the absolute values of the smallest Dirac eigenvalues are bounded from below at the scale . The current mass should be selected somewhat above this scale, in order to cause the chiral symmetry breaking. The reason is obvious from the expression for the quark condensate, as the trace of the massive quark propagator
(37) 
Tuning the quark mass values used in the simulation is important: it should be as small as possible, yet large enough to avoid large finitevolume effects.
An example of how this tuning works in practice is demonstrated in Fig.4. In the left panel we show the Dirac eigenvalue distributions for 6 values of the current quark mass. As one can see, the spectrum itself changes, from a dip near at zero mass, to a sharp peak at larger masses (as for quenched case). The CasherBanks “density of the eigenvalues at the origin” thus changes dramatically, from zero to very large values: hence, reading off is clearly not a good method for estimating the quark condensate estimate. By contrast, the points shown in the right panel of Fig.4 show the quark condensate estimated from the loop expression (37). We observe that the quark condensate remains approximately mass independent for larger masses , while decreasing for smaller masses, eventually going to zero at .
Let us first comment on the continuous lines in Fig.4(b), passing through our 6 simulation points. Those correspond to the loop expression (37) in which one uses the eigenvalues as determined from the simulation with the corresponding dynamical mass, but in the loop formula itself the mass is substituted by a separate parameter, called the “valence quark mass” in literature. Such “valence approximation” has been used before, sometimes for quenched ensembles, to approximate the mass dependence outside of the region where simulations were done. However, as seen from the plot, the trend indicated by those curves is quite misleading, being completely different from what is obtained from the dynamical simulations – the points themselves. It shows that the main mass dependence is stored in the eigenvalue distribution, reflecting the dynamical quark mass of the vacuum, of the “sea” quarks.
Let us now remind the SmilgaStern theorem Smilga:1993in , which is based on chiral perturbation theory and predicts the slope of the eigenvalue spectrum near the origin, as a function of the chiral condensate , the pion decay constant and the number of flavors:
(38)  
Lattice and instanton liquid simulations had indeed confirmed that the slope of the distribution changes sign at , from a negative constant at to positive at .
Let us also point out that the spectral density in the expanded form near can be put into the integral in (37), generating the condensate expanded in powers of . The zeroth order is the CasherBanks term. Returning to Fig.4(b) for , we see that no liner term in is expected, and higher order terms can be neglected: so the correct infinite volume limit should be just a constant. Our guess of its value is indicated by a horizontal dashed line. Comparing it to the data we see that for finite solum effects are large and those should not be used. Since we cannot afford to do simulations with variable masses in all cases studied, we will below only use the dynamical mass , at which we expect the condensate accuracy determination to be not worse than say 20%.
Our computed eigenvalue spectra, for different volumes and number of flavors, are shown in Fig.5. While some finitevolume effect is still seen at the smallest eigenvalues, it is much less prominent than in zero mass simulations discussed above: they are only seen for .
The spectra are clearly “gapped”, without any small eigenvalues, and the gap magnitude grows as the dyon density falls. We thus conclude that case needs rather large dyon density to break the chiral symmetry, perhaps close to the highest one we used.
Other simulation show no gap and thus correspond to chirally broken phase, except for very low dyon density. At different dyon density and number of flavor, we performed a linear fit . The fit was restricted to a range with until the maximum for which the reduced remains smaller than 1. This way, we identified constants and with the quark condensate and the SmilgaStern constant, respectively. The result of the fit is shown in Fig.6. The nonzero condensate observed for the and ensembles at high dyon density is rather densityindependent.
More systematic account for the finitevolume effects and more accurate determination of the condensate can be done by comparing our data with the expectations for the so called “mesoscopic regime”, in which the volume is not macroscopically large. Quantitative predictions of the shape are known from chiral random matrix models Shuryak:1992pi , which are well confirmed on the lattice. We plan to do so elsewhere.
v.4 Free energy, dyons’ interaction, the backreaction and confinement
The well known perturbative holonomy potential, which has a minimum at zero holonomy, has been argued by Polyakov’s original work to get cancelled by the nonperturbative effects, resulting in a vanishing Polyakov line and thus confinement. Lattice study of the effective holonomy potential, in particularly recent Ref.Diakonov:2012dx , have indeed found such behavior. Diakonov Diakonov:2009ln further argued that the nonperturbative contribution comes from the back reaction of the dyons. He also argued that (at least the Coulomblike moduli part of) the dyon interaction is small and can be approximately ignored. In this section we make the first step toward understanding of whether those ideas are correct.
We calculate the change in the free energy of the system, between the interacting and noninteracting dyon gas using standard a thermodynamic integration method, based on the adiabatic switching of the interaction between the particles. One splits the action into independent particles and their interaction, introducing the adiabatic parameter
(39) 
and do simulations with changing from 0 to 1. The resulting freeenergy is recovered as follows
(40) 
The resulting dependence of the action as a function of the adiabatic parameter is shown in Fig.7. Of course, we include all three ingredients of the interaction, the moduli, screening and the fermions, and measure them separately as well: but for brevity we will not discuss those details here.
As one simulates an ensemble with a fixed number of dyons, one can disregard any constant factors in the partition function, as the relative weights of the ensemble configurations depends only on terms which are functions of the dyonic collective coordinates. This explains why one has positive freeenergy of the noninteracting ensemble at , which is not really physical and depends of what factors we do or do not include in the free gas expression. The freeenergy change is physical, and a shift to negative is typical for liquids. The value itself is for the entire system of 64 dyons: it thus corresponds to a free energy change per particle of which is due to the attraction in the system.
To put it in perspective, one should know the action per dyon. We will discuss available lattice data in the Appendix: the effective action fitted to those are about per dyon. The result of this section show that the average interaction part of it is reasonably small, . On the other hand, the partition function is enhanced by about per dyon, not a negligible enhancement.
As for the back reaction to the total holonomy potential, we think that it is premature to discuss it now, before lattice simulations provide more define estimates on the dyon density.
Vi Summary and further developments
This paper is the first study of the topology, in the finitetemperature QCD near , which incorporates the nonzero holonomy and the (anti)self dual dyons as the ingredients. We (i) formulate the partition function for (anti) selfdual dyons incorporating oneloop screening, moduli space measure and sermonic zero modes; and (ii) perform the first direct simulation of the resulting statistical ensembles. We obtain certain set of results, showing how the dyon system depends on the dyon density (changed by about 2 orders of magnitude), on the number of light quark flavors , from 0 to 4, and on the (current) quark mass . We have in particularly studied the Dirac eigenvalue spectra and chiral properties: the gaps in the symmetric phase and the quark condensate in the broken phase. More generally, we may say that in spite of the longrange nature of Coulomb forces involved and linearly rising screening, the simulations of the dyonic system turned out to be not as challenging as we thought when we started this work. While the correlations between constituents are locally strong, overall liquidlike behavior is observed.
At the moment we are not calculating the fermionic propagators and correlation functions of operators with various quantum numbers, but it can be done rather directly, in a way similar to that in the instanton liquid.
Mapping those results to finiteT QCD can be done provided the information about the dyons, such as their density we reviews briefly in the Appendix, are extracted from the lattice. One should also directly relate fermionic quasi zero modes with the underlying topology, which we believe will display the dyonic dominance in the chiral symmetry breaking phenomenon, in the near region. One may also expect progress in respect with the back reaction of the dyons on the holonomy potential and the origin of confinement.
Another direction is to extend the model toward theories in which one has highT confinement under the analytical control, see e.g. Poppitz:2012nz , by including adjoint fermions with modified (e.g. periodic or with an arbitrary phase) boundary conditions on the torus.
Extending our results to lower temperatures – or much smaller holonomy values – remains a numerically challenging task.
Attractive correlations of the LM pairs should mostly recover the original instantons, yet preserving
small Polyakov line value and confinement.
At the same time, fermions induce the attractive interaction, producing topologically neutral clusters
whose role in the QCD vacuum has not been elucidated or studied.
Acknowledgements The first acknowledgements on this subject obviously goes to Pierre van Baal and his students, who pioneered the instantontodyon transition at nonzero holonomy. Another source of inspiration came from beautiful works of Dmitry Diakonov, who sadly is no longer with us. ES acknowledges the role of his collaborator and former student Tin Sulejmanpasic. The work of ES is supported by the USDOE grant DEFG88ER40388. The work of PF was partially performed at the Interdisciplinary Laboratory for Computational Science (LISC), a joint venture of Trento University and Bruno Kessler Fundation.
Appendix A Distances on the sphere
In our program, all distances are defined on a 3dimensional sphere embedded in 4 dimensions. Let be a unitary vector in 4 dimension, which will be associated to the position of the pseudoparticle of type in the th molecule. Monte Carlo trial moves are made by performing random rotations of given small angle (boldness) on the sphere.
The angle between the 4D vector identifying the pseudoparticle of type relative to the molecule and that points to the pseudoparticle of type relative to the molecule on the unitary sphere is given by:
(41) 
(Notice that this formula is illdefined for . On the other hand, the distance of a particle from itself is never needed in the calculation.) Hence the distance between these pseudo particle is simply:
(42) 
where is a length scale which can be identified with the box size.
For illustration purposes we mapped the to a flat space using the stereographic projection. The three dimensional coordinates of our dyons are defined on a sphere, hence we use 4dim vectors with the constraint .
One may transform them into the 3 polar angle angles like this
(43) 
The inverse formulas are:
(44) 
In a stereographic projection to the flat 3dspace in polar coordinates the and angle keep the same meaning as above, while the distance from the origin is defined by the angle according to
(45) 
(R=1 is the sphere radius and the argument of tan is the angle conjugated to psi in the equilateral triangle of the projection).
Note that maps the south pole to the origin , and the equator is mapped to . Hence the cartesian coordinates on the projected flat space are:
(46) 
Appendix B The Coulomb fields on a
The usual coordinates have the following line element (metric)
(47) 
Standard Poisson equation is
(48) 
For a point charge at the north pole () the solution is obviously independent on and we need to solve only the dependent part
(49) 
For zero r.h.s. one get solution which looks like a positive charge at one pole and the negative charge at the other. Indeed, on a compact manyfold all field lines have to go somewhere, and symmetry require it to be the opposite pole.
One remedy can be introducing the homogeneous compensating charge. This solution is
(50) 
It is compared with the naive Coulomblike and solution given above in Fig. One can see from it that its behavior is very reasonable, and it vanishes at the opposite pole.
Screening effect contains the averaging of the dipole potential over the 3d space, which produces the linear rising effective potential (as explained in the text). Keeping one charge at the north pole and spreading another at 2d sphere one gets another solution of the Poisson eqn. It is
(51) 
with . Integrating the square of this field over the 3d sphere one gets the following expression
(52) 
Both agree reasonably well till about the equator of the 3d sphere.
Appendix C SU(2) holonomy on the lattice
One important potential input is the holonomy value (i.e. the VEV of the Polyakov line’s VEV ). While we have not yet used it, we indicate that its temperature dependence for the SU(2) gauge group is known. To our knowledge the most accurately determined values are from Ref. Hubner:2008ef . The resulting fit of the Polyakov line VEV
(53) 
where is the nearcritical parameter
(54) 
The temperature dependence of the holonomy is then .
Note that these data and the fit are consistent with the indices of the Ising universality class to which the SU(2) pure gauge theory belongs.
Appendix D The (anti)selfdual dyons on the lattice
While we only consider a case with 2 colors, we remind the reader that one of the major benefits of the dyon approach is the smooth behavior of the action in the large limit. Indeed, the dyon partition function is roughly the root of the instanton partition function, so the dyon action in the large limit is expressed via ’t Hooft coupling as
(55) 
While at very strong coupling it is small, at the “physical” value of the parameter usually used for strongly coupled QGP applications it is still in the range of 34 units. Hence, the weight factor is still a rather small number, providing a substantial suppression of the dyon density. It is important to stress that, although in principle it may be argued that the value of the action is not large enough to justify the usage of the semiclassical approach, the corresponding dyon ensemble remains dilute even close to , hence providing an a posteriori justification of the approach.
While we generally postpone mapping of the results of our dyonic model to the QCD thermodynamics directly, in this section we still discuss some available results on the dyons, from Refs. Ilgenfritz:2006ju and Bornyakov:2008im . One of the issues here is whether they are or are not compatible with the with semiclassical expressions, in spite of possible low accuracy of those.
These works both investigate the SU(2) gluodynamics on the and lattices respectively, adopting the same treelevel improved Symanzik action at a bare coupling . In the former lattice, this choice corresponds to the critical temperature, , while in the latter lattice to a temperature .
Note that lattice practitioners often introduce adhoc units of energy and length which are defined in such a way to mimic the dynamics of the physical world. For example, In order to fix an energy scale, they measure the dimensionless constant.
(56) 
and artificially set the vacuum string tension value to , i.e. to its physical value in the ’realworld’ QCD. Furthermore, one keeps also the ’realworld’ relation between the length and the energy units . Under this set of conventions, one obtains the following “absolute” values for the lattice spacing and the critical temperature
(57)  
(58) 
When comparing to lattice results from those works – such as the density of various dyons – we may adopt such absolute units. However, in general, we will refrain from doing so because this may lead to confusion when studying the dependence on the number of flavors . In addition, holding fixed the value of the string tension constant is only practical in a lattice approach, in which this observable can be straightforwardly calculated. In view of such considerations, in the present analysis we avoid using “physical units” and consider only dimensionless ratios and couplings.
In Ref.Ilgenfritz:2006ju , Ilgenfritz and coworkers report a list of the absolute densities of static (type) dyons and calorons ( dyon pairs) which were found in their lattice calculation. As dyons are much lighter than , the calorons density is dominated by the density of the heaviest dyons. Their results are plotted in Fig.8. In a more recent paper by Bornyakov and coworkers Bornyakov:2008im an improved identification of the dyon type was preformed by adopting the fermion method and also by computing the value of the Polyakov loop at the dyon center: it should be 1 for the dyons, but for “additionally twisted” dyons.
Their result is shown in the same figure by the filled diamond, and its suggests the density about twice larger than in the other method. Unfortunately, the latest study was done only at one temperature . The factor two discrepancy between the results in Ref. Ilgenfritz:2006ju and in Ref. Bornyakov:2008im is due to different efficiency of their topologysearching algorithms. The values of the the paper Ref.Ilgenfritz:2006ju mean to be the lower bound on the densities, as the authors also report also rather large number of “unidentified topological clusters” in the same table.
It is instructive to compare the dyon densities observed on the lattice with the analytic semi classical expressions. Those include the holonomy and the effective action in a combination
(59) 
where is the normalization constant and the product of the preexponent powers of the holonomy is taken from the DGPS weight for the LM pair Diakonov:2004jn . The density of dyons is given by obvious change .Since the holonomy goes to small values at high more rapidly than the charge is running, one expects a significant excess of the dyons over the type ones at high . The difference between the two densities should however disappear at , with a characteristic singularity related to critical behavior of the SU(2) theory.
As it is clear from Fig.8 , the observed dependence shown by curves can be well reproduced by the semiclassical expression, with two fitted parameters:
(60) 
We also note that the point from Bornyakov:2008im at corresponds to rather small dimensionlless density
(61) 
or volumes/dyon above even the highest values used in our simulations.
We recall that the lattice data apply to a fully interacting ensemble, e.g. including the screening part of the action and that the data/expressions discussed so far apply only to a pure gauge theory, i.e. to , which is characterized by the weakest coupling along the critical line. Hence, by increasing number of quark flavors we expect to assess a significantly stronger coupling domain, hence a smaller . If so, eventually there would be no numerical suppression and the diluteness should become . Such a dense liquid of (nonsemiclassical) dyons is obviously not expected to be described by the simple formulas used above. This regime should be the subject of future theoretical and numerical studies.
While grouping of the dyons into instantons suggests about equal number of dyons, this equality is only reached at and below , while at the lighter dyons are more numerous than heavier ones. The relevent lattice data are shown in Fig.8. As one can see, the density of dyons grows at high temperature, as their action decreases. The density of type, associated with calorons, decreases, at large . We emphasize again that the identified topological objects in the lattice simulations are just a subset of the of the observed topological clusters.
Only the solid diamond corresponding to the analysis of Bornyakov et al. Bornyakov:2008im should be considered as a quantitative determination of the type dyon density. Indeed, this was done by combining the fermionic quasizeromodes filtering procedure with the action smearing to locate topology, and by computing the sign of the Polyakov line at the center to identify the and type dyons. Such a work focused on SU(2) gauge configurations in a lattice. These authors confirmed that several main predictions of the semiclassical theory remain valid in their lattice configurations. In particular, all fermionic zero and nearzeromodes are locally chiral ( a general property of their topological nature), the periodic fermions interact with and antiperiodic fermions with dyons, as expected from semiclassical solutions. The dyons are dilute and paired with ones into clusters” discussed in Shuryak:2012aa . The chiral symmetry of the antiperiodic fermions is unbroken. The type dyon ensemble is dense and thus chiral symmetry for periodic fermions is broken. They also found that while the Ltype dyons occupy only about 3% of the lattice volume, the M dyons occupy a significant fraction of it.
References
 (1) S. Mandelstam, Phys. Rept. 23, 245 (1976).
 (2) G. ’t Hooft, Nucl. Phys. B 190, 455 (1981).
 (3) C. Gattringer, Phys. Rev. Lett. 88, 221601 (2002).
 (4) E. V. Shuryak, Nucl. Phys. B 203, 93 (1982).
 (5) D. Diakonov and V. Y. .Petrov, Nucl. Phys. B 272, 457 (1986).
 (6) L. Y. .Glozman, C. B. Lang and M. Schrock, arXiv:1207.7323 [heplat].
 (7) T. Schafer and E. V. Shuryak, Rev. Mod. Phys. 70, 323 (1998) [arXiv:hepph/9610451].
 (8) P. Faccioli, A. Schwenk and E.V. Shuryak, Phys. Rev. D. 67, 11 (2003).
 (9) M. Cristoforetti, P. Faccioli and M. Traini, Phys. Rev. D 75, 054024 (2007) [hepph/0701223 [HEPPH]].
 (10) M. Cristoforetti, P. Faccioli, M. C. Traini and J. W. Negele, Phys. Rev. D 75, 034008 (2007) [hepph/0605256].
 (11) P. Faccioli, Phys. Rev. D 65, 094014 (2002) [hepph/0111322].
 (12) M. Cristoforetti, P. Faccioli, G. Ripka and M. Traini, Phys. Rev. D 71, 114010 (2005) [hepph/0410304].
 (13) P. Faccioli, A. Schwenk and E. V. Shuryak, Phys. Lett. B 549, 93 (2002) [hepph/0205307].
 (14) P. Faccioli and E. V. Shuryak, Phys. Rev. D 65, 076002 (2002) [hepph/0108062].
 (15) P. Faccioli, Phys. Rev. C 69, 065211 (2004) [hepph/0312019].
 (16) M. Cristoforetti, P. Faccioli, E. V. Shuryak and M. Traini, Phys. Rev. D 70, 054016 (2004) [hepph/0402180].
 (17) N. I. Kochelev and V. Vento, Phys. Rev. Lett. 87, 111601 (2001).
 (18) K. M. Lee and P. Yi, Phys. Rev. D 56, 3711 (1997) [arXiv:hepth/9702107].
 (19) T. C. Kraan and P. van Baal, Phys. Lett. B 435, 389 (1998) [arXiv:hepth/9806034].
 (20) F. Bruckmann, D. Nogradi and P. van Baal, Acta Phys. Polon. B 34, 5717 (2003) [arXiv:hepth/0309008].
 (21) D. Diakonov, Topology and Confinement, arxiv:0906.2456v1 (2009)
 (22) E. Shuryak and T. Sulejmanpasic, Phys. Rev. D 86, 036001 (2012) [arXiv:1201.5624 [hepph]].
 (23) E. Shuryak and T. Sulejmanpasic, The Chiral Symmetry Breaking/Restoration in Dyonic Vacuum II, The adjoint fermions. in progress
 (24) D. Diakonov, N. Gromov, V. Petrov and S. Slizovskiy, Phys. Rev. D 70, 036003 (2004) [arXiv:hepth/0404042].
 (25) M.F.Atiyah and N.J.Hitchin, The Geometry and Dynamics of Magnetic Monopoles, Princeton University press 1988.
 (26) G.W.Gibbons and N.S. Manton, Nucl. Phys. B274 (1986) 183.
 (27) E. V. Shuryak, Sov. Phys. JETP 47, 212219 (1978)
 (28) R. D. Pisarski, L. G. Yaffe, Phys. Lett. B97, 110112 (1980).
 (29) E. M. Ilgenfritz, B. V. Martemyanov, M. MullerPreussker and A. I. Veselov, Phys. Rev. D 73, 094509 (2006) [heplat/0602002].
 (30) V. G. Bornyakov, E. M. Ilgenfritz, B. V. Martemyanov and M. MullerPreussker, Phys. Rev. D 79, 034506 (2009) [arXiv:0809.2142 [heplat]].
 (31) A. V. Smilga and J. Stern, Phys. Lett. B 318, 531 (1993).
 (32) K. Huebner and C. Pica, PoS LATTICE 2008, 197 (2008) [arXiv:0809.3933 [heplat]].
 (33) E. V. Shuryak and J. J. M. Verbaarschot, Nucl. Phys. A 560, 306 (1993) [hepth/9212088].
 (34) D. Diakonov, C. Gattringer and H. P. Schadler, JHEP 1208, 128 (2012) [arXiv:1205.4768 [heplat]].
 (35) E. Poppitz, T. Schaefer and M. Unsal, arXiv:1212.1238 [hepth].