Deconfinement Phase Transition in the SU(3) Instanton-dyon Ensemble
DDeconfinement Phase Transition in the SU (3) Instanton-dyon Ensemble
Dallas DeMartini and Edward Shuryak
Center for Nuclear Theory, Department of Physics and Astronomy,Stony Brook University, Stony Brook NY 11794-3800, USA
Confinement remains one the most interesting and challenging nonperturbative phenomenon innon-Abelian gauge theories. Recent semiclassical (for SU(2)) and lattice (for QCD) studies havesuggested that confinement arises from interactions of statistical ensembles of instanton-dyons withthe Polyakov loop. In this work, we extend studies of semiclassical ensemble of dyons to the SU (3)Yang-Mills theory. We find that such interactions do generate the expected first-order deconfinementphase transition. The properties of the ensemble, including correlations and topological susceptibil-ity, are studied over a range of temperatures above and below T c . Additionally, the dyon ensemble isstudied in the Yang-Mills theory containing an extra trace-deformation term. It is shown that sucha term can cause the theory to remain confined and even retain the same topological observables athigh temperatures. I. INTRODUCTION
Quantum Chromodynamics (QCD) is the quantumfield theory describing the fundamental particles andforces that make up nuclear physics. While QCD is re-markably successful in describing nuclear physics, manyphenomena remain beyond the scope of what can be stud-ied analytically. Notably, nonperturbative phenomenasuch as confinement – the disappearance of quarks andgluons from the physical spectrum – is not completely un-derstood. Confinement occurs not just in QCD, but invarious Yang-Mills theories with or without quarks, mak-ing it clear that it emerges from the non-perturbative be-havior of the gluons, rather than the quarks. Above cer-tain critical temperature T c deconfinement takes place,and the QCD-like theories turn into a new form of mat-ter, the Quark-Gluon Plasma (QGP).Historically the first mechanism of the deconfinementtransition was a ’dual superconductor’ model [1–3]. At T < T c the chromoelectrically-charged quarks and glu-ons are connected by QCD flux tubes, dual to magneticflux tubes in superconductor. With the advent of lat-tice gauge theories many aspects of this scenario wereput to the test. In particular, the profile of the QCDflux tubes [4] was found to agree well with dual super-conductor model. Monopoles were observed and foundto rotate around these flux tubes, as expected. Bose-Einstein condensation of monopoles was detected and itscritical temperature was shown to coincide with T c [5].A high density of monopoles was found to be responsiblefor unusual kinetic properties of QGP [6].Euclidean formulation of the gauge theory lead to dis-covery of 4D topological solitons known as BPST instan-tons [7]. A model of their ensemble, the Instanton LiquidModel (ILM) [8], has explained how instantons generatechiral symmetry breaking. As certain extrema of thepath integral over gauge configurations, they form a ba-sis for semiclassical theory, consistently including fluctua-tions around classical fields. Furthermore, one can studythe interaction between instantons in their statistical en-semble: those studies explained behavior of correlation functions of various mesonic and baryonic currents, forreview see e.g. Ref. [9]. Yet the instanton theory has notreproduced confinement.Euclidean formulation of finite temperature QCD nat-urally led to a nonzero value of the Polyakov loop (cid:104) P (cid:105) (cid:54) = 0as a signature of deconfinement. Note that through-out this paper we use (cid:104) P (cid:105) as shorthand for (cid:104) T r [ P ( (cid:126)x )] (cid:105) .It can be interpreted as the nonzero vacuum expecta-tion value (VEV) of the time component of the gaugefield A , also known as nonzero holonomy . The natu-ral question was then how to deform the instanton con-figurations in a way consistent with nonzero holonomyin the bulk. It was answered in Refs. [10, 11], whodiscovered that instantons dissolve into N c (number ofcolors) constituent solitons, called instanton-dyons (orinstanton-monopoles). Like original instantons, they are(anti)selfdual, so their actions and topological charges areequal. But, unlike instantons, their actions and topologi-cal charges are not quantized to integers; standard indextheorems are avoided because instanton-dyons have mag-netic charges and therefore are still connected by Diracstrings.It was then realized that instanton-dyons provide avery valuable bridge between the theory of monopolesand instantons, providing a way to explain both confine-ment and chiral symmetry breaking in a single setting.Unlike monopoles, the instanton-dyons are semiclassicalobjects, allowing for the construction of a consistent the-ory of an interacting ensemble. Last but not least, thestatistical sum in terms of instanton-dyons is ”Poissondual” to that based on monopoles, see Refs. [12, 13].Semiclassical approaches to finite- T gauge theories,with and without quarks, have lately been subject ofmultiple studies. Confinement in this theory is due tothe back reaction of the dyon ensemble on the Polyakovloop, forcing it to take zero value at T < T c , see e.g. Ref.[14] for a simple model, Ref. [15] for mean-field analysis,and Refs. [16, 17] for numerical simulations of the SU (2)gauge theory. It is important that they are semiclas-sical objects, unlike the QCD monopoles [13], with theactions S dyons ∼ /g ∼ log ( T /
Λ) growing with temper-ature. Therefore their densities are suppressed at high T a r X i v : . [ h e p - ph ] F e b as an inverse power of T . At temperatures comparable tothe critical one, T ∼ T c , numerically S dyons / (cid:126) ≈
4. Theinverse of it is the small parameter of our semiclassicalexpansion.Although in this work we study pure gauge theory,rather than QCD-like theories with light quarks, let usmention that the instanton-dyon ensemble also describesthe breaking of chiral symmetry. Numerical studies ofchiral phase transition can be found in Refs. [18, 19].Using light quark Dirac eigenstates, one can identify theindividual dyons inside lattice configurations from large-scale QCD simulations. As shown recently [20], the low-est Dirac states are indeed generated by the instanton-dyons, in the form remarkably insensitive to large densityof perturbative quarks and gluons in which they are im-mersed.In this work we extend the studies [16, 17] to the pure SU (3) gauge theory. Instead of two species of instanton-dyons, we now deal with three. It is well-known fromlattice studies [21] that the pure SU (3) gauge theory(and all pure SU ( N c ) theories with N c ≥
3) possessesa first-order phase transition, rather than the second-order transition seen in SU (2). So the question we ad-dress below is whether this, and other, features of pure SU (3) gauge theory can or cannot be reproduced by theinstanton-dyon ensembles.This paper is structured as follows: The dyon interac-tions and partition function are discussed in Section II.Section III lays out some technical details of the simula-tion as well as the data analysis performed. The physicalresults, including correlation functions and the tempera-ture dependence of parameters such as the Polyakov loopare shown in Section IV. Finally, in Section V, the de-confinement transition of the dyon ensemble is studied inthe trace-deformed Yang-Mills theory. II. INSTANTON-DYONS IN THE SU (3) GAUGETHEORYA. Dyons and holonomy
The Polyakov loop is an order parameter of the decon-finement phase transition P = P exp( i (cid:73) A aµ T a dx µ ) , (1)where T a is the color generator in the fundamental rep-resentation. Its VEV, (cid:104) P (cid:105) is a unitary matrix withphases for eigenvalues, so A can be defined as A =2 πT diag ( µ , µ , µ ). These phases are related to theholonomies of the dyons by ν i = µ i +1 − µ i , where µ i +1 > µ i and µ = µ . The connection between thePolyakov loop and confinement can be seen through itsrelation to the free energy of a static quark (cid:104) P (cid:105) = e − F q /T . (2) Clearly, (cid:104) P (cid:105) = 0 is the confining holonomy, correspond-ing to F q → ∞ , removing massive quarks from the spec-trum completely. By combining the previous definitionsof the holonomies and phases, the relationship betweenthe holonomy and the average Polyakov line is seen to be (cid:104) P (cid:105) = 13 + 23 cos(2 πν ) . (3)From this it is clear that ν = 1 / (cid:104) P (cid:105) =0),while ν → (cid:104) P (cid:105) = 1).For SU ( N c ) gauge theories, there are N c types ofdyons: N c − M i -type dyons, theycorrespond to the maximal diagonal subgroup. One moretype is called the ”twisted” (that is, 4-time-dependent) or L -type dyon. There are corresponding antidyons as wellgiving in total 2 N c species of dyons. Note that we will use species to refer to one of the six dyons and type to refer toone of the three pairs of a dyon and its antidyon (i.e. an M -type dyon refers to both M and ¯ M dyon species).The action of each individual dyon of type i is denotedby S ν i , where ν i are the so-called holonomy parameters.Those are also called “fractions of the holonomy circle”because they satisfy the sum rule N c (cid:88) i ν i = 1and generically have N c − r i ∼ / (2 πν i T ).In SU ( N c ) theories, (anti)instantons are comprised ofone of each of the types of (anti)dyons, although thisis not necessarily the cause for other gauge groups. Forexample, in SU (3), the instanton I = M LM and theantiinstanton ¯ I = ¯ M ¯ L ¯ M .The main parameter defining the properties of thedyons is the classical instanton action S = 8 π g ∼ log ( T )containing the coupling at the temperature-dependentscale. While this parameter is independent of theholonomy, the actions of individual dyons of type i are S i = S ν i . In the SU (2) case, one defines one holon-omy parameter ν ∈ [0 , M dyon having action S M = S ν and the L dyon having the conjugate action S L = S (1 − ν ). In the SU (3) case there are two diago-nal Gell-Mann matrices, λ and λ and in principle twoholonomy parameters. However, the average Polyakovloop (cid:104) P (cid:105) is a gauge-invariant, and thus physical, object.Restricting it to be real enforces ν M = ν M reducingthe holonomy to just a single parameter for SU (3) with ν ∈ [0 , / S M = S M = S ν, S L = S (1 − ν ) . (4)This symmetry between the M - and M -type dyonsmeans that they have equal densities as well, greatly re-ducing the space of parameters that needs to be consid-ered. µ µ µ ⌫ ⌫⌫ L M M FIG. 1. (Color online) Structure of the holonomies and dyon types in SU (3). Circle has circumference of 1. there is a perturbative interaction between thermal gluons and the holonomy. This generates the Gross-Pisarski-Ya↵epotential V GP Y [22], which appears in the an exponent in the partition function. For SU (3), this is given by V GP Y T ˜ V = 4 ⇡ ⌫ (1 ⌫ )) + (2 ⌫ (1 ⌫ )) ) . (5)Here the potential is shown with factors divided out so that it has units of free energy density. This potential disfavorsconfinement, having minimum at trivial holonomy ⌫ = 0, and a maximum at the confining holonomy ⌫ = .The dyon contribution to the partition function consists of two parts, the contributions of the dyons in the absenceof interactions Z , which can be expressed analytically directly from the input parameters, and Z np the contributionsof the dyon interaction, whose calculation is the subject of the simulation of the partition function, performed byMonte-Carlo algorithms in this work. The total partition function is their product Z = Z GP Y Z Z np .The statistical weight for a single instanton was explicitly calculated in Ref. [23] for SU (2). Taking the dilute limit,which removes interaction between the dyons, this gives Z SU (2)0 = ⇤(4 ⇡ ) S e S ⌫ ⌫ (1 ⌫ ) ⌫ )3 (6)It is easy to see how this factors into the weight for each individual dyon: Each dyon species contributes a factor of p ⇤ S e p S o ⌫ and the two holonomy terms stem from the holonomies of the individual M- and L-type dyons. A factorof 4 ⇡⌫ i is divided out for each dyon type. This is to remove a constant term that appears remains when taking thedilute limit in Z np . Knowing this, it is easy to construct the instanton weight in SU (3) and extend it to an arbitrarynumber of dyons by summation. The partition function, assuming an equal number of dyons and antidyons, is Z = X N M ,N L ,N M ✓ N M ! ( ˜ V d ⌫ ) N M ◆ ⇥ ✓ N L ! ( ˜ V d ⌫ ) N L ◆ ⇥ ✓ N M ! ( ˜ V d ⌫ ) N M ◆ , (7)where d ⌫ is the weight of an individual dyon with holonomy ⌫ , d ⌫ = ⇤4 ⇡ S e S ⌫ ⌫ ⌫ . (8)Now we may take the limit V ! 1 and assume equal densities of M - and M -type dyons. The free energy FIG. 1. (Color online) Structure of the holonomies and dyontypes in SU (3). Circle has circumference of 1. B. The partition function and dyon interactions
A complete calculation of the dyons’ free energy re-quires the construction of the dyonic partition function.We start first with effects that are not induced by thedyons’ non-perturbative interactions. In the absence ofall dyonic effects, there is a perturbative interaction be-tween thermal gluons and the holonomy. This generatesthe Gross-Pisarski-Yaffe potential V GP Y [22], which ap-pears in the an exponent in the partition function. For SU (3), this is given by V GP Y T ˜ V = 4 π ν (1 − ν )) + (2 ν (1 − ν )) ) . (5)Here the potential is shown with factors divided out sothat it has units of free energy density. This potentialdisfavors confinement, having a minimum at the trivialholonomy ν = 0, and a maximum at the confining holon-omy ν = .The dyon contribution to the partition function con-sists of two parts, the contributions of the dyons in theabsence of interactions Z , which can be expressed ana-lytically directly from the input parameters, and Z np thecontributions of the dyon interaction, whose calculationis the subject of the simulation of the partition function,performed by Monte-Carlo algorithms in this work. Thetotal partition function is their product Z = Z Z np .The statistical weight for a single instanton was explic-itly calculated in Ref. [23] for SU (2). Taking the dilutelimit, which removes interaction between the dyons, thisgives Z SU (2)0 = Λ(4 π ) S e − S ν ν − (1 − ν ) − ν )3 − (6)It is easy to see how this factors into the weight for eachindividual dyon: Each dyon species contributes a factor of √ Λ S e −√ S o ν and the two holonomy terms stem fromthe holonomies of the individual M- and L-type dyons.A factor of 4 πν i is divided out for each dyon type. Thisis to remove a constant term that appears remains whentaking the dilute limit in Z np . Knowing this, it is easyto construct the instanton weight in SU (3) and extendit to an arbitrary number of dyons by summation. Thepartition function, assuming an equal number of dyonsand antidyons, is Z = (cid:88) N M ,N L ,N M (cid:18) N M ! ( ˜ V d ν ) N M (cid:19) × (cid:18) N L ! ( ˜ V d − ν ) N L (cid:19) × (cid:18) N M ! ( ˜ V d ν ) N M (cid:19) , (7)where d ν is the weight of an individual dyon with holon-omy ν , d ν = Λ4 π S e − S ν ν ν − . (8)Now we may take the limit V → ∞ and assume equaldensities of M - and M -type dyons. Additionally weresolve the factorial terms with Stirling’s approxima-tion carried out to three terms ln N ! ≈ N ln N − N +ln( √ πN ). The free energy F = − ln( Z ) is given by thefollowing expression f = 4 π ν (1 − ν )) + (2 ν (1 − ν )) ) − n M ln (cid:20) d ν en M (cid:21) − n L ln (cid:20) d − ν en L (cid:21) + ln(8 π N M N L )˜ V + ∆ f (9)where ∆ f is the free energy density stemming from theinteractions of the dyons. If the dyons have classical bi-nary interactions ∆ S class and a volume metric G , theircontributions to the partition function and the free en-ergy density are Z np = 1˜ V (4 N M +2 N L )3 (cid:90) Dx det ( G ) e − ∆ S class (10)∆ f = − ln( Z np ) (11)The set of parameters that minimizes the free energydensity corresponds to the physical dyon ensemble in theinfinite volume limit. This elucidates the main procedureof this work: to first compute the free energy density ofthe ensemble for a wide range of parameters, and thenlocate their values which minimize it.The classical interactions between dyons and antidyonsare, at distances exceeding the dyon cores, asymptot-ically Coulomb-like. For generic SU ( N c ) theories, theinteractions are given in Ref. [24]. At shorter lengthscales, the interaction is modified. The interaction be-tween dyons and antidyons of the same type was studiedin detail and parameterized by Larsen and Shuryak [25].In the previous SU (2) model, this interaction was usedfor dyon-antidyon pairs of the same type, while dyonsand antidyons of different types had purely Coulomb-likeinteractions. In this work we use a single modified pa-rameterization for all dyon-antidyon pairs. The parame-terization used is given by∆ S d ¯ dclass = − S C d ¯ d π ( 1 rT − . π √ ν i ν j e − . π √ ν i ν j rT ) , (12)where C d ¯ d is a coefficient with value 2 for pairs of thesame type and − SU (2), two mod-ifications have been made: the substitution ν i → √ ν i ν j was made to accommodate pairs of dyons with differentholonomies (e.g. M ¯ L ), and the coefficient in front ofthe exponential term has been reduced from the originalvalues of 3 . r = x / (2 πνT ). Dyons pairs ofthe same type (regardless of duality) experience a repul-sive core. While the core potential has not been studiedin detail it reasonably described by∆ S coreclass = νV e πνT ( r − r ) . (13)The two parameters in this interaction, x and V arechosen on a phenomenological basis and should be sub-ject to constraints from appropriate lattice data whenpossible. See Appendix B for a discussion of their valuesand the effects of changing them.Additionally, the dyons experience an effective poten-tial from the fluctuation determinant of the instanton.This effect leads to the Diakonov determinant [26] of themetric of the space of dyons’ collective variables G im,jn = δ ij δ mn (4 πν m − (cid:88) k (cid:54) = i T | r i,m − r k,m | + (cid:88) k T | r i,m − r k,p (cid:54) = m | )+ 2 δ mn T | r i,m − r j,n | − δ m (cid:54) = n T | r i,m − r j,n | , (14)where r i,m is the position of the i’th dyon of type m.This metric only accounts for the selfdual dyons. Anequivalent metric ¯ G is used for the antidyons as well.The metric is true for dyons of different species at anydistance, but only true at larger distances for dyons ofthe same species. We modify the terms with a cutoff r → (cid:112) r + (3 / πT ) so that the diagonal entries go to0 rather than −∞ when ν = . The effective potentialfrom this metric comes as∆ S − loop = − ln(det G det ¯ G ) . (15) The diagonal of this metric does not vanish in the dilutelimit, instead going to (cid:81) i πν i . This has been accountedfor by modifying the individual dyon weight in Eq. (8).It is observed that in the case the densities of differentdyon types are unequal, then in the infinite volume limit,the sum diverges. We therefore regulate all Coulombterms in both the classical and one-loop potentials withthe dimensionless Debye mass r → re M D rT . Like thecore parameters, the value of the Debye mass is a choicethat should be subject to improvement.Now let us give a brief qualitative discussion of howsuch interactions should generate confinement. The maininteraction driving confinement is the repulsive cores ofthe dyons. Because the volume of the cores goes as1 /ν i , this interaction disfavors any type of dyon becom-ing large and thus favors the confined phase ν = whereall dyons are the same size. The long-distance classi-cal Coulomb interactions disfavor confinement however.Because dyons repel anti-dyons of different types but at-tract antidyons of the same type, these interactions fa-vor ensembles with many dyons of the same kind - alarge value of n M /n L - rather than equal numbers of allspecies. At larger values of n M /n L the entropy portionof the partition function Z more strongly favors smallvalues of ν , driving the system to the deconfined phase.The one-loop Coulomb-like interactions have the oppositesign and thus opposite effect. These one-loop interactionsare suppressed by a factor of S and are only comparableto the classical interactions at low T .The nonperturbative interactions grow weak comparedto the perturbative portions of the partition function as T increases. At low temperature we have an ensem-ble where the cores and one-loop interactions dominate,binding the dyons into instantons and at high tempera-tures we have a gas of mostly M i ¯ M i pairs bound by theirclassical attraction and the entropy and Gross-Pisarski-Yaffe terms have shifted ν to a lower value. III. THE SIMULATION SETTING AND DATAANALYSIS
Numerical integration over the dyons’ coordinates isperformed by Monte-Carlo techniques. The position ofeach dyon is updated and then accepted or rejected ac-cording to the standard Metropolis algorithm. Once thisis done for each dyon five times (five updates is approx-imately the autocorrelation time of the ensemble), theconfiguration is sampled and its properties are computed.The free energy is obtained via the usual integration overauxiliary parameter λ : e − F ( λ ) /T = (cid:90) Dre − λS , (16)∆ F = (cid:90) dλT (cid:104) ∆ S (cid:105) , (17)where this nonperturbative free energy is added to theperturbative contribution from Eq. 9. Integration overthe dummy parameter is performed in 10 equal steps λ = 0 .
1, ... 1, with the average energy at each stepbeing computed from 2000 configurations. Because thecontribution at small λ is large, the value of the aver-age interaction at λ = 0 is obtained from a linear fit ofthe nearest points and included in the numerical integra-tion. With all of this, the statistical uncertainty in thecomputations of f are generally at the percent level.The goal of these simulations is to calculate the freeenergy density f of ensembles with three main input pa-rameters, the holonomy ν and the dyon densities n M , n L over a range of temperatures near T c . All simulationsare performed with a fixed number of dyons N D = 120or 122 (60 or 61 dyons and 60 or 61 antidyons). For thesix species of dyons in SU (3), N D = 4 N M + 2 N L . Thesystem is studied on flat geometry with finite density im-posed by a periodic box surrounded by 26 image boxes,as was done for SU (2) in Ref. [17]. The densities of thedyons, n M and n L are controlled by varying the size ofthe box and the ratio of the numbers of each type of dyon N M /N L for a specific value of N D . Two values of N D areused to allow for more values of N M /N L to be used.There are two technical details of the simulation thatmust be mentioned. The first is the Diakonov determi-nant det G , which is not positive definite for randomly-placed dyons. It has been shown by Bruckmann et. al.[27] that in fact, at reasonable densities, nearly all ran-dom configurations have at least one negative eigenvalue λ i . This is handled by a modification of the effectivepotential∆ S − loop = (cid:26) − ln(det G det ¯ G ) if λ i ≥ ∀ iV max else (18)where V max is a large value (we choose V max = 100)which serves to reject all Metropolis updates into theregion of negative eigenvalues. The choice of V max isarbitrary so long at it is large enough to ensure thatno configuration with negative eigenvalues is acceptedduring the course of the simulation.The second is an approximation which must be madedue to computational resources. The total system, themain box and all its nearest images, contains 27 N D dyons. Its 2D version is shown in Fig. 2. The mostdemanding step in the simulation is the computation ofdet G which is an O ( N ) process. Performing this com-putation at every update step of the Metropolis algo-rithm is numerically costly, given the large number ofparameters for which the free energy must be calculated.Instead, the change in action associated with each up-date is computed considering only the dyons in a box ofthe same size as the simulation box centered at the po-sition of the dyon being updated. This approximation iswell-justified by the fact that all long-range interactionsare regulated by a Debye screening length.Once the free energy density of each set of parametersis determined, a fit is performed to determine the depen- FIG. 2. (Color online) Schematic depiction of the main simulation box (center) and surrounding images setup. Dashed boxshows the volume containing all dyons whose interactions are considered in the Metropolis update step of the dyon representedby a (red dashed) square. Such a volume always contains N D dyons. in the simulation is the computation of det G which is an O ( N ) process. Performing this computation at everyupdate step of the Metropolis algorithm is numerically costly, given the large number of parameters for which the freeenergy must be calculated. Instead, the change in action associated with each update is computed considering onlythe dyons in a box of the same size as the simulation box centered at the position of the dyon being updated. Thisapproximation is well-justified by the fact that all long-range interactions are regulated by a Debye screening length.Once the free energy density of each set of parameters is determined, a fit is performed to determine the dependenceof the free energy density near the minimum. For each value of the action S there is a 3-dimensional space ofinput parameters with coordinates ( ⌫, n M , n L ). The points near a minimum are fit to a quadratic function in thesecoordinates. From the fit, the minimum value of f and its location in the coordinate space is determined. For everyvalue of S tested, all local minima are identified and the points near them are fit as described. The fit with the lowestminimum value of f is the global minimum. For cases where a minimum is located in the in the confined phase with ⌫ = and n M = n L , these equalities are assumed exact and the minimum is fit to a quadratic function only in n M .In summary, the results presented in this work are obtained by the following procedure:(i) perform simulations over a wide range of parameters ⌫ , n M , n L for each temperature tested,(ii) determine the global minimum for each temperature via fitting to determine the physical values of these parameters ⌫ ( T ), n M ( T ), n L ( T ),(iii) study the ensembles with the physical values with increased statistics to compute correlations and topologicalsusceptibility.(iv) With these increased-statistics runs, we generate 50,000[?] configurations with the physical interactions betweenthem ( = 1, we no longer need to integrate over it) using the fitted values of the inputs for each temperature studied.In order to more accurately reproduce the fitted values of both dyon densities simultaneously, the number of dyons isallowed to vary more with 126 N D FIG. 2. (Color online) Schematic depiction of the main sim-ulation box (center) and surrounding images setup. Dashedbox shows the volume containing all dyons whose interactionsare considered in the Metropolis update step of the dyon rep-resented by a square point. Such a volume always contains N D dyons. dence of the free energy density near the minimum. Foreach value of the action S there is a 3-dimensional spaceof input parameters with coordinates ( ν, n M , n L ). Thepoints near a minimum are fit to a quadratic functionin these coordinates. From the fit, the minimum valueof f and its location in the coordinate space is deter-mined. For every value of S tested, all local minima areidentified and the points near them are fit as described.The fit with the lowest minimum value of f is the globalminimum. For cases where a minimum is located in thein the confined phase with ν = and n M = n L , theseequalities are assumed exact and the minimum is fit to aquadratic function only in n M . TABLE I. Input parameters used for the main simulationruns. Some additions and changes were made as needed.(Some non-integer values of S were added near T c and thevalues of n M were increased at lower temperatures.)min. max. step size no. of steps S ν . n M N M /N L In summary, the results presented in this work are ob-tained by the following procedure:(i) perform simulations over a wide range of parameters ν , n M , n L for each temperature tested, FIG. 3. (Color online) Holonomy dependence of the free energy density for different values of the M -dyon density n M in bothphases. Left: confined phase with S = 12 . n L = n M Right: deconfined phase with S = 13 . n L = n M / . ¯6. Error barsnot shown for readability. (ii) determine the global minimum for each temperaturevia fitting to determine the physical values of these pa-rameters ν ( T ), n M ( T ), n L ( T ),(iii) study the ensembles with the physical values withincreased statistics to compute correlations and topolog- ical charge distributions.(iv) With these increased-statistics runs, we generate80,000 configurations with the physical interactions be-tween them ( λ = 1, we no longer need to integrate overit) using the fitted values of the inputs for each temper-ature studied. IV. RESULTSA. Structure of the holonomy potential and thephase transition
It is expected that the pure SU (3) gauge theory pos-sesses a first-order phase transition at some critical tem-perature T c . The dyons interactions generate the holon-omy potential which has a global minimum in the con-fined phase at T < T c and then a jump to a new globalminimum in the deconfined phase at T > T c . The firstgoal of this work was to check for the existence of such aphase transition. Parameters of the interaction were cho-sen so this occurs around S ∼
12 (See Appendix B formore details). For each value of the temperature, which isrelated to S by Eq. (A1), the shape of the holonomy po-tential was computed and the minimum was determined.The value of the holonomy and dyon densities at whichthis minimum occurs are the physical values the ensembletakes on at that temperature.The instanton action S was varied in steps of 1 from S = 8 - 21. It was observed that the phase transitionoccurs between 13 < S <
14 and additional simulationswere performed at S = 12 .
75, 13 .
25, 13 .
5, and 13 .
75 tolocate the critical temperature more accurately. A linearfit was performed to the free energy density with the nearest points on both sides of the transition to determinethat the phase boundary is located at S = 13 .
18. Withthis, the action can now be rewritten in terms of
T /T c .The range of temperatures studied here is then 0 . T c 530 MeV.This model of the dyon interactions relies on a semi-classical expansion. An important consideration is then,whether or not the main semiclassical quantity, the actionof the dyons S i = S ν i , remains sufficiently large overthis temperature range. In the confined phase, all dyonshave action S / 3, meaning there is clearly a lower boundwhere the action becomes small and we choose S = 8, S i = 2 . ¯6 as a reasonable lower bound. In the decon-fined phase, as the instanton action grows, the holonomydecreases making the action of the M i -type dyons thelimiting factor. The action of these dyons remains nearlyconstant at S M (cid:39) T c . Thus, the dyons remainsufficiently semiclassical throughout the entire range oftemperatures studied.By performing similar linear fits, all parameters of theensemble can be determined on both sides of the phasetransition. Table II summarizes these values. From thistable it is clear that the transition is first order - theholonomy / Polyakov loop and both dyon densities arediscontinuous at the phase boundary. TABLE II. Values of the parameters of the ensemble aboveand below the critical temperature T c from linear fits to thenearest data points on either side of the phase transition at S = 13 . T → T − c T → T + c ν / (cid:104) P (cid:105) n M n L At temperatures below T c the ensemble is in the con-fined phase with ν = 1 / n M = n L , as can be seenin Fig. 3 (left). At densities below the physical one,the minimum shifts to the left as the nonperturbativeinteractions become weak compared to the perturbativecontribution to the free energy. At higher densities theensemble prefers to remain in the confined phase, butwith a larger free energy minimum and curvature of thepotential.The deconfined phase has a similar structure, but withthe global minimum occurring at ν < / n M > n L .However at densities above that of the global minimum,the minima continue to move towards larger ν . At thesedensities, the repulsive core dominates and it becomes en-ergetically favorable to make the many M i dyons smallerat the cost of making the few L dyons larger. It is possi-ble that at densities higher than what were studied here,the ensemble may have a minimum at ν > / S for specific valuesof n M /n L . The structure of the first-order phase transi-tion can be seen more clearly in Fig. 4. By consideringthe minimum free energy density selected from all com-binations of n M and n L as a function of the holonomy,the two local minima – one in the confined phase andthe other in the deconfined – are clearly visible. At thevalue of S nearest to the critical value, the free energyat the two minima are nearly degenerate and the globalminimum switches between the two as the temperaturechanges.This structure is different from the that of SU (2). In SU (2), where the phase transition is second order, ratherthan having two degenerate minima, the holonomy po-tential flattens near T c (see e.g. Fig. 5 of Ref. [16]). Thisallows the minimum to quickly, but smoothly, shift fromthe confining holonomy to smaller values. Additionally,there is a ν ↔ − ν symmetry not present in SU (3). B. Temperature dependence of the parameters The free energy density f , unlike other physical quan-tities, remains continuous across even a first-order phasetransition, as we see in Fig. 5. Its derivative, however,may not. The free energy varies with temperature muchmore rapidly in the confined phase than the deconfined. FIG. 4. (Color online) Holonomy dependence of the minimumfree energy density near the phase transition. Error bars notshown for readability.FIG. 5. Temperature dependence of the free energy densityof the dyon ensemble. The most important feature of the dyon ensemblefor describing the deconfinement transition is the av-erage Polyakov loop as a function of the temperature (cid:104) P ( T ) (cid:105) . Below T c , the holonomy takes the confiningvalue ν = 1 / (cid:104) P (cid:105) = 0. At T c the value jumps to ∼ . T increases. The valueof the average Polyakov loop above the phase transitionshows qualitative agreement with the lattice data [21],but does not increase with temperature as quickly. A FIG. 6. (Color online) Temperature dependence of the aver-age Polyakov loop of the dyon ensemble. Lattice data takenfrom Ref. [21] and shown without error bars. Error on latticedata has magnitudes comparable to dyon data. change to the parameters of the dyon interactions couldimprove the agreement with the lattice data.Another set of properties of the ensemble are the den-sities of each dyon type, shown in Fig. 7. At T < T c ,all densities are equal reflecting the fact that, in the con-fined phase, all dyons have equal statistical weights, coresizes, and masses and have symmetric interactions be-tween them. At T > T c , the holonomy decreases discon-tinuously, causing a similar change in both dyon densi-ties. In the case of the M -dyon density n M , the decreasein ν increases the statistical weight of the M i -type dyons,while simultaneously increasing the size of their repulsivecores. These competing effects result in what is seen inthe ensemble: a small, O (5%), decrease in n M at thephase transition. The L -dyon density sees a more sig-nificant decrease due to the large decrease in statisticalweight of the L dyons from the increase in 1 − ν .Overall, the dyon densities decrease with temperature,consistent with the expected behavior of instantons. Itis the dyon densities that set the upper limit on the tem-perature that can be studied in our ensemble. As tem-perature increases, the ratio n M /n L becomes larger andlarger, and, at the largest values, the number of L dyonsin the simulations is just 1 or 2. The small number of L dyons makes fitting near the free energy minimum in n L difficult as their contribution to f is vanishingly small.Thus, our upper limit ( S = 21) is a technical constraint,rather than a physical one; accurately probing higher S would require a larger ensemble with better statistics. FIG. 7. (Color online) Temperature dependence of the densi-ties of each type of dyon in the ensemble. C. Spatial correlations between dyons It is useful to study the effects of the dyons’ interac-tions by studying the spatial correlations between them.Such correlations will be useful for comparison to stud-ies of the dyons in lattice configurations. The moststraightforward characteristic of the spatial correlationsare the correlation functions C ij ( rT ) between two dyonsof species i and j . In SU (3), because the instanton iscomprised of three constituent dyons, it is useful to alsodefine the ’instanton correlation function’ C I ( yT ), where y is the hyperdistance in the 6-dimensional space of Ja-cobi coordinates of the three dyons y = 13 (cid:0) r M ,L + r M ,M + r L,M (cid:1) . (19)In both cases, the functions are normalized such that C ij , C I → rT ) for the two-dyon cor-relation functions and ( yT ) for the instanton correlationfunction).Many of the correlation functions are redundant so weneed to only observe a subset of the correlation func-tions to understand the behavior of the ensemble. Fig.8 shows these functions in both phases. The strongestcorrelation between the dyons is seen in the instantonchannel where the three constituent dyons have positiveshort-range correlation due to the attractive terms in theDiakonov determinant. FIG. 8. (Color online) Spatial correlations of the dyons in the M ¯ M (upper right), L ¯ L (upper middle), M ¯ L (upper right), M M (lower left), LL (lower middle), and instanton (lower right) channels above and below the phase transition. In several channels we can see clearly the effects ofthe repulsive core at short ranges. In the D ¯ D channels,we can see positive correlation just beyond the cores.Changes in the cores can also be seen: in the M ¯ M channel, the core becomes larger and softer above T c ,soft enough that many pairs can still be found at x < x ,while the cores of the L dyons becomes smaller andharder, reducing the correlation function to 0 at x < x .The smaller core size also allows the attractive long-rangeinteractions to become larger near to the core, enhancingthe correlation found there.Despite having no long-range interactions, there issmall correlation seen in the same-species ( M M and LL ) channels beyond the core due to some mutual corre-lations to other dyons. Compared the the D ¯ D channelsthe anticorrelation at small r is even stronger due to ad-ditional repulsion from the Diakonov determinant. Somecorrelation functions for the dyons in pure SU (3) Yang-Mills theory at T /T c = 0 . 96 have been observed in latticeconfigurations [28]. We find qualitative agreement withthese results in all channels except for dyons of the samespecies. Here some short-range correlation is observedamong the lattice dyons, suggesting an absence of therepulsive core used in our model in this channel. D. The vacuum angle and moments of thetopological charge Another set of important characteristics of the dyonensemble are its topological properties. Non-abeliangauge configurations generically may possess some topo- logical charge Q with its density q ( x ) given by q ( x ) = g π (cid:15) µνρσ F aµν ( x ) F aρσ ( x ) . (20)The (anti)instantons are topologically nontrivial objectswith Q = ± 1. The constituent dyons then each carrysome fraction of this charge. Dyons of species i carrytopological charge Q = ν i while antidyons have charge Q = − ν i . Because we only consider ensembles withequal numbers of dyons and antidyons, (cid:104) Q (cid:105) = 0 in anysub-volume of the box. The most prominent topologicalfeature of the gauge theory, which has been extensivelystudied on the lattice [29–31], is the topological suscep-tibility χ = (cid:104) Q (cid:105) /V .Because the dyons in our semiclassical ensemble havewell-defined positions and topological charges, measur-ing the average topological charge only requires summingover the charges of the dyons in a given sub-volume. Thisis done by splitting the main simulation box in half alongeach of the 3 axes and computing the total charge inthe said half-boxes. This results in 3 independent mea-surements per configuration, meaning the average of anypower of Q is computed from 240,000 measurements.The temperature dependence of the topological sus-ceptibility over the entire temperature range probed inthis work is shown by blue points in Fig. 10. To makecomparison to lattice data easier, we show the relative susceptibility, normalized all sets to their values just be-low the critical temperature. The issue of absolute com-parison is further discussed in Appendix C.With χ ( T ), like (cid:104) P ( T ) (cid:105) , we see a clear first-order phasetransition in the data. Unlike with the Polyakov loophowever, such an abrupt transition with a jump is not0 FIG. 9. Distribution of values of the total topological charge Q in the half-boxes at S = 12 ( T /T c = 0 . obvious from the lattice data points (although they ofcourse do not contradict existence of a jump). The topo-logical observables have been analyzed on the lattice inmany works, we compare our dyon results to those inRefs. [29, 32], which cover a similar temperature rangeas we have studied. FIG. 10. (Color online) Topological susceptibility as functionof temperature relative to the topological susceptibility justbelow T c . Lattice data taken from Refs. [29, 32]. Latticevalues are normalized relative to their stated T = 0 values.Data from Ref. [29] (green triangles) is so-called ’2-smear’data. Error bars on dyon data are smaller than the pointsize. The (Euclidean) Lagrangians of generic SU ( N ) gaugetheories can be appended by an additional topologicalterm L = 14 F aµν ( x ) F aµν ( x ) − iθq ( x ) , (21)where θ is the so-called vacuum angle and q ( x ) is thetopological charge density as previously defined in Eq.(20). A non-zero θ explicitly breaks CP symmetry. InQCD it is known that | θ QCD | (cid:47) − , as it is con-strained by the upper bound on measurements of theneutron’s electric dipole moment [33, 34].By expanding the free energy density f ( θ ) around θ =0, its dependence on θ can be studied at small, but non-zero θ . It can be expanded as [35] f ( θ ) = f (0) + 12 χθ (1 + b θ + b θ + ... ) , (22)where the coefficients b n are related to cumulants of thetopological charge computed at θ = 0. We compute thefirst two terms, which are given explicitly by b = − (cid:104) Q (cid:105) − (cid:104) Q (cid:105) (cid:104) Q (cid:105) b = (cid:104) Q (cid:105) − (cid:104) Q (cid:105)(cid:104) Q (cid:105) + 30 (cid:104) Q (cid:105) (cid:104) Q (cid:105) . (23)In the confined phase, we find that all b values areconsistent with b being constant below T c . These sevenvalues vary around 0 . − . 03 and have an average b ( T < T c ) = 0 . T c , b quickly drops andthen remains approximately constant just above 0 . T = 0 value of b has been particularly wellstudied on the lattice [36–38] with values around b (0) = − . b approaches the value predicted by theDilute Instanton Gas Approximation (DIGA), b ( T ) →− / T c , our values of b are consistent with the mag-nitude of the T = 0 value predicted on the lattice, butwith opposite sign. In the high temperature limit, weagain see that our dyon model predicts positive valuesrather than negative, as well as being an order of magni-tude smaller than the lattice results. Clearly, these smallnon-Gaussianities in the topological charge distributionare quite sensitive to the dyon interactions. Changes tothe dyon interactions could improve the agreement withlattice data.Finally, our values for b are compatible with zero forall temperatures as is also observed on the lattice [36, 40].Upper limits from these lattice results constrain its valueto be | b ( T = 0) | < − .1 FIG. 11. Expansion coefficient b as a function of tempera-ture. V. DYONS IN THE TRACE-DEFORMEDGAUGE THEORYA. Polyakov loop and the phase transition Trace-deformed gauge theories (TDGT) were intro-duced in Refs. [24, 41], adding certain terms containingpowers of the Polyakov line to the theory action, such as∆ S def = h (cid:90) d x | P ( (cid:126)x ) | (24)with a new parameter h . At T > T c , when in the un-deformed theory the Polyakov loop is nonzero, the newterm obtains an additional contribution and suppressesits value, shifting the theory back toward the confiningholonomy. For large enough h the system returns to the (cid:104) P (cid:105) = 0 phase, which we will call the “reconfined phase”.These theories have been studied in numerous latticesimulations, of which we will mention those of the Pisagroup [42–44]. Their main findings are that for a specifictemperature above T c at large enough h > h ∗ ∼ / 3, atwhich (cid:104) P (cid:105) returns to zero, there is no more dependenceon h , and the reconfined phase is remarkably similar tothe original (low temperature) confined phase. It wasshown to possess the same spectrum, topological observ-ables, and even spectrum of periodic strings (torelons).From the perspective of the instanton-dyon theory thatwe are developing, it is clear that the reconfined phase,like the confined one, corresponds to symmetric valueof the holonomy ν = 1 / L, M , M of the dyons. Yet the temperaturescale, and thus the periodicity interval of the Euclidean time τ , is different. Therefore, the absolute scale of thedyon action parameter S = 8 π /g ( T ) must increase,making the system more dilute. And yet, the topolog-ical susceptibility χ is the same as at T = 0 despite asignificantly reduced dyon density!It is possible to study the instanton-dyon ensemble in atheory with arbitrary h by the addition of a deformationterm to the free energy density ∆ f def = h (cid:104) P (cid:105) . Becausethis new term depends only on the value of the holonomy,one does not need to perform new MC simulations tostudy the effect of changing h ; one needs only to add thenew term to the un-deformed results and perform newfits to find the minima.The first questions to ask are how does changing h affect the value of the Polyakov loop and at what valueof h does the system become reconfined? Fig. 12 showsthe temperature-dependence of the Polyakov loop for afew different values of h . As expected, the larger thetrace deformation, the more suppressed the values of (cid:104) P (cid:105) become at T > T c . Additionally, we see that increasing h increases the critical temperature of the theory as well,as the Polyakov loop is suppressed back to the confiningvalue. FIG. 12. (Color online) The value of the Polyakov loop asa function of temperature (cid:104) P ( T ) (cid:105) for different values of thetrace-deformation parameter h . The critical temperature of the trace-deformed theoriesare determined via the same method as the original the-ory: the intersection of linear fits to f on both sides pf thetransition. At this point we are limited to determining T c for h ≤ 40, as above this value the the Polyakov loop re-mains very close to the confining value and the minimumvalue of the holonomy fit become closer to ν = thanthe holonomy step size used for the simulations, meaning2that ν = is compatible with our uncertainties for alltemperatures studied. In order to accurately determinethe critical temperature for large h , smaller steps in both S and ν are necessary. This also makes determining thenature of the phase transition inconclusive at larger h .At small h , it is clear to see that the phase transitionremains first order, while at larger h (above ∼ 10) we donot have sufficient resolution in ν to determine whetherthe holonomy potential continues to have two distinctminima or becomes a smooth crossover. FIG. 13. Critical temperature of the theory with trace-deformation parameter h relative to that of the the pure Yang-Mills theory. While we are currently limited in how high of a criticaltemperature we can probe (up to ∼ . T c (0) as per Fig.13), the theory may have non-trivial behavior at highertemperatures. In the limit that T → ∞ the density ofdyons goes to zero and only two terms in the free en-ergy density remain: the (deconfinement-favoring) Gross-Pisarski-Yaffe potential and the (confinement-favoring)trace-deformation term. Clearly which term dominatesdepends on the value of h . At h > π / (cid:39) . 74 theconfining holonomy is the global minimum.This leads to the conclusion that there are three dis-tinct regimes of the TDGT:i) h < π / 18: The theory is confined at low tempera-tures and deconfined at high temperatures.ii) 5 π / < h < h max : The theory is confined at lowtemperatures, deconfined in some intermediate regionand then confines again at high temperatures.iii) h > h max : The theory is confined at all temperatures.While we can’t yet access high enough temperaturesto directly see the second phase transition in regime (ii),the last value of (cid:104) P (cid:105) for h = 40 in Fig. 12 suggests that (cid:104) P (cid:105) may be starting to decrease and return to zero. B. Topological observables of the reconfined theory One of the most interesting features of the deformedgauge theory observed on the lattice [42–44] is that whenthe trace deformation is large enough to return the sys-tem to the confining holonomy, the topological observ-ables, namely χ and b , also return to the same valuesthey had in the confined phase of the un-deformed theory,and then show no more dependence on further increasing h . Rather than choose specific temperatures and vary h aswas done in the lattice studies, we choose an arbitrarily-large value of h and vary temperature. As mentionedin the previous section, for some large value of h , suchthat h > h max , the system should remain in the con-fined phase at all temperatures. In our dyon model, wecan study this scenario, which we call the ’maximally-deformed’ theory, by simply demanding that ν = anddetermining the dyon densities that minimize the free en-ergy density with this constraint. Because our individualsimulations have fixed holonomy value ν = , (cid:104) P (cid:105) = 0the trace-deformation term does not contribute to thefree energy density of these simulations and doesn’t needto be considered in the fits. FIG. 14. Dyon density of the maximally-deformed theoryas a function of temperature. Both dyon densities are equal n = n M = n L . The question then is whether or not the dyon ensem-ble, like the lattice theory, possesses the same topologicalobservables in the maximally-deformed theory as in theoriginal confined phase. In the case of χ , our dyon modeldid not predict a constant value but an approximatelylinearly-decreasing one. We observe analogous behaviorto what is seen on the lattice; in the maximally-deformed3 FIG. 15. Topological susceptibility of the maximally-deformed theory as a function of temperature. Error barsare smaller than the point size. theory the abrupt jump in χ disappears and it displaysthe same temperature dependence above T c (0) as below.Comparing the dyon density (Fig. 14) and the topologi-cal susceptibility (Fig. 15), one can see that the decreasein χ is mostly driven by the decreasing dyon density. Al-though n i decreases by a factor of ∼ χ decreasesby a factor of ∼ (cid:104) Q (cid:105) decreases by a factor of ∼ b is more inconclusive. It remainsconstant up to about 1 . T c (0) and then decreases quicklyand is compatible with zero at all temperatures abovethat point. Again, it remains positive at all tempera-tures. It should be noted that the error bars on the b measurements are still quite large and the trend couldlook quite different with improved statistics.In order to determine what drives this behavior, it isuseful to look at some correlation functions in this de-formed theory. It is clear that the ensemble has somenon-trivial temperature dependence even though it re-mains in the confined phase. The two main channels withpositive correlation are, just like the un-deformed the-ory, the D ¯ D channels and the instanton ( M M L ) chan-nel. As the temperature is increased, the ensemble morestrongly prefers D ¯ D pairs rather than (anti)instanton ascan be seen in Fig. 17. This is due to the increase inthe strength of the D ¯ D classical attraction. Eventually,these pair correlations grow strong enough to completelykill correlation in the instanton channel.In the original theory, the system prefers binding intoinstantons in the confined phase and then transitions into FIG. 16. Expansion parameter b in the maximally-deformedtheory as a function of temperature. preferring M ¯ M - and M ¯ M pairs at high temperatures;the density of L ¯ L is becoming small at this point. How-ever, because the core sizes of the M i dyons grow large,the binding in these channels does not become too large.The maximally-deformed theory then displays behav-ior not seen in the original theory. At high temperaturethe system becomes a three-component gas of D ¯ D pairs,with all types having equal densities and core sizes. Be-cause the theory remains confined, the core sizes neverbecome large and the correlation between dyons in D ¯ D pairs becomes very large. The topological observables arevery sensitive to this binding , as the (anti)instantons and D ¯ D pairs have topological charges Q = ± Q = 0,respectively. The shift towards the latter leads to a neu-tralization of the topological charge and is responsible fordriving the values of χ and b down at higher tempera-tures.Looking to the future, one can see that this behaviorhas interesting implications when quarks are added tothe theory. It is known that chiral symmetry restorationis directly related to the quark zero modes on the in-stantons. In particular, the quark interactions drive thesystem to form instanton-antiinstanton ’molecules’. Thisneutralizes the fluctuations of the topological charge andshifts the Dirac eigenvalues away from zero at high tem-peratures, restoring chiral symmetry [45]. In fact in thiswork we already observe D ¯ D pairing at higher temper-atures. We see here that such correlations continue togrow, even when the theory remains confined due to thedeformation term in the action. This suggests that inthe trace-deformed theory with quarks one may see chiralsymmetry restoration occur inside the confined phase.4 FIG. 17. (Color online) Spatial correlation functions of the dyons in the M ¯ M (left) and instanton (right) channels in themaximally-deformed theory at three different temperatures. VI. SUMMARY AND DISCUSSION This paper reports the first direct numerical simu-lations of the ensemble of instanton-dyons in the pure SU (3) Yang-Mills theory. We use classical and semi-classical one-loop measures which have derived the inter-dyon interactions. We have performed Monte-Carlo sim-ulations of an ensemble of dyons in a 3D cube with peri-odic boundary conditions over a range of temperatures,densities, and holonomy values.Our main objective is to see whether the semiclassicalensemble of instanton-dyons does or does not reproducethe deconfinement transition, which is, in this case, ofthe first order. We indeed find that there are two distinctphases, the confined (with free energy minimal at holon-omy value 1 / 3, Fig. 3 (left)) and the deconfined (withfree energy minimal elsewhere, Fig. 3 (right)). We wereable to map out the holonomy potential and see explic-itly the two nearly-degenerate minima of the potentialrepresenting the transition between the two phases (redsquares in Fig. 4) and study various properties near thecritical temperature.Comparison between our semiclassical model and lat-tice data for VEV of the Polyakov loop (cid:104) P ( T ) (cid:105) is shownin Fig. 6. Remarkably, the jump of it at T c is reproducedrather precisely, with some moderate deviations at highertemperatures. As shown in Fig. 7, this large jump resultsin L -dyon suppression in the deconfined phase, while thedensity of M i -dyons has a rather smooth temperaturedependence.The pair-wise spatial correlations between variousdyons are shown in Fig. 8. Overall, they are in qual-itative agreement with expectations. Note also, that for a system with a strong first order transition, one does notsee drastic changes in the correlations, except for overallchange of scale due to the change in the temperature.We have studied fluctuations of the topological charge.Since our simulations each have a fixed number of dyons,this is done by counting charges in the half-box. A typ-ical distribution is shown in Fig. 9, at first sight havingjust a normal Gaussian form. Its width - the topologi-cal susceptibility χ - is compared to lattice data in Fig.10. The magnitude of the jump and behavior in the de-confinement phase are found to be quite similar, but thedependence at T < T c is somewhat different. Our resultsfor non-Gaussianity parameter b are perhaps still toonoisy. Let us note that these observables, as well as thedyon correlations, are very sensitive to the details of thedyon interactions, particularly the short-range classicalinteractions, which are the strongest. These interactionsare at present the most poorly-understood aspect of thedyon model with some parts being analytically derived(e.g. the Diakonov determinant) and others being simplyphenomenological (e.g. the repulsive core). In order toachieve a quantitative agreement between the dyons andlattice data, if possible, these interactions should be thesubject of further study and improvement.We also report the first study of trace-deformed gaugetheory in the instanton-dyon setting. As expected, byincreasing the coupling h of the deformation term, oneindeed suppresses the jump and value of the Polyakovloop in the deconfined phase. The phase transition loca-tion is also pushed to higher values, see Fig. 13. Overall,like on the lattice, we see that increasing h leads to areturn to the topological susceptibility as in the confin-ing phase of the un-deformed theory. We however still5see that this agreement is only partial, see for examplethe spatial correlations taken in the maximally-deformedconfining theory.Completing this summary of our results, we again em-phasize that the semiclassical model used, which hasmany orders of magnitude fewer degrees of freedom thanlattice gauge theory, is able to explain properties of thedeconfinement phase transition of the SU (3) gauge the-ory quite well. Some quantities – like the jump in thePolyakov line – are reproduced precisely, others semi-qualitatively, but we have not seen any serious disagree-ments. The double-trace deformation study also showsthat both our model and lattice react to this addition tothe action in a very similar way.Thinking of the future, the next step will of course bethe inclusion of dynamical quark flavors via their zeromodes and interactions. Such work will allow not onlyfor the study of confinement (which will now be expectedto be a smooth crossover), but also to study the role ofthe dyons in chiral symmetry restoration. ACKNOWLEDGMENTS This work is supported by the Office of Science, U.S.Department of Energy under Contract No. DE-FG-88ER40388. The authors also thank the Stony BrookInstitute for Advanced Computational Science for pro-viding computer time on the SeaWulf computing cluster. Appendix A: Units Temperature is the main physical quantity of thiswork: it defines the holonomy via A = 2 πνT as wellas the physical sizes of the dyons. As seen in previ-ous equations, all interactions of the dyons contain thedimensionless distance rT rather than the dimensionfuldistance r . Because of this we define many quantities interms of this dimensionless distance: the dimensionless3-volume ˜ V = V T , the dyon densities n i = N i ˜ V , andthe free energy density f = F ˜ V T .All parameters are presented dependent on the instan-ton action S which is related to the temperature by S = 8 π g = ( 113 N c − N f ) ln (cid:18) T Λ (cid:19) (A1)where, in pure SU (3) gauge theory, N c = 3 and N f = 0.As with previous work on the ensemble of SU (2) dyons[16], the value of the parameter Λ is a choice that mapsthe instanton action S to the temperature T . This pa-rameter was chosen to be Λ = 2 . 8. It should be thesubject of improvement in future comparison with lat-tice data. The reader should keep in mind that this valueis different than the one chosen for the previous SU (2)studies, Λ SU (2) = 1 . 5. Because the temperature scaleis set by this choice , conversion to physical units may be done setting the dyons’ critical temperature equal tothat of the lattice gauge theory: T c = 9 . ↔ 260 MeV[46, 47]. Appendix B: Comments on interaction parametersand finite-size effects1. Parameters of the dyon interactions Let us take a moment to remind the reader thatthis is a model of the dyon ensemble. Certain param-eters of the dyon interaction, namely V , x , and M D ,are not known from first principles, and are thus phe-nomenological parameters of the model. In the previ-ous SU (2) work [17] these parameters were chosen to be( V , x , M D ) = (20 , . , . SU (3) model, it was found that these parameterswere unsuitable for generating the desired phase transi-tion. These parameters are sensitive to the details of thetheory and one should not expect that they are identicalfor all N c .In this work, the parameters used were ( V , x , M D ) =(10 , / , . S ∼ 12, jumping from (cid:104) P (cid:105) = 0to (cid:104) P (cid:105) (cid:39) . 4. The choice of x was not arbitrary, as whenthe factor of 1 /ν is included, the dyons in the confinedphase have the same dimensionless core radius as wasused in SU (2). This choice of parameters exists on a 3D’island’ in the parameters space in which such choices ofthe parameters possess reasonable properties. They canand should be further constrained by future data regard-ing identification of the dyons in lattice configurations. Inparticular, if dyon densities and correlations can be stud-ied with sufficient statistics, such data can reveal detailsof the dyons’ interactions. Some correlations betweendyons on the lattice (in physical QCD) can be found inRef. [48], although the statistics are not yet sufficient tomake such constraints.Each of these parameters has a clear effect on the phasetransition of the ensemble. Large values of any of theseparameters favors the confined phase and drives T c up-ward. In the case of V and x , large values make therepulsive cores larger, more strongly favoring the confin-ing holonomy. A large value of M D suppresses the large-range interactions making the dyon cores comparativelymore important. For a more systematic analysis of theeffects of changing the parameters see Ref. [17].Additionally, the Debye mass is not necessarily aconstant, but rather a temperature-dependent quantity M D ( T ). In fact, in the previous work by one of the au-thors [16] the Debye mass was treated as an input param-eter to be swept over, much like the densities and onlyconfigurations with input values consistent with the itsdefinition M D = g V ∂ F∂ν | n , (B1)6were chosen. This increases the number of configurationsto be tested, but allows for a more general, temperature-dependent Debye mass. Both methods have been shownto produce reasonable properties for the SU (2) dyon en-semble, so we have saved computational time by treatingit solely as an input. One could also consider a more’direct’ method in the future, simply taking the value of M D ( T ) as an input from lattice data. 2. Finite-size effects In principle, equilibrium statistical mechanics is de-rived in the thermodynamic limit N D , ˜ V → ∞ , which isnot possible in numerical simulations. For such a finitenumerical simulation, the main concern is whether thesystem size used is sufficiently large enough so that themeasured observables can be reliably extrapolated to thethermodynamical limit. The simplest way to check theseeffects is to simulate larger systems with the same pa-rameters and study the dependence of some observableson the size N D .Here we discuss a representative example of the finite-size effects as seen in Fig. 18. In increasing the size ofthe system by up to a factor of 2, there is a noticeablechange in the free energy landscape. Looking at Eq. (9),there are only two terms which can vary with system size:the dyon interaction term ∆ f and the third term stem-ming from Stirling’s approximation, which has explicit˜ V -dependence. Increasing the system size decreases thecontribution from the third term as it goes as 1 / ˜ V . Sub-tracting this term from the results reveals that the contri-bution from the dyon interactions increases with systemsize. Both of these effects are reduced in the deconfinedphase where the density is lower.The main factor in the finite-size effects regarding dyoninteractions comes from dyons near the boundaries of thesystem. As long as a finite number of dyons and imageboxes are used, there are dyons near the faces of thecube whose total short-range interactions are reduced.The fraction of dyons near the faces of the cube shouldscale as N − / D and thus vanishes in the thermodynamiclimit. Parameters that determine the effective range ofinteractions such as the core sizes and Debye mass alsoplay an important role in determining how quickly theresults converge to the N D → ∞ limit.Different quantities may converge at different speeds.From Fig. 18 it is seen that the value of the free energyminimum decreases by ∼ 2% when doubling the numberof dyons. A simple linear fit of f as a function of 1 /N D gives an infinite-volume value of f = − . . 5% de-crease from the N D = 120 value. The input parametersof the ensemble - the holonomy and densities - dependonly on the location of the minimum, which clearly variesmuch less. Fitting the curves to quadratic forms givesthe densities n M = 0 . . . 607 for dyon numbers N D = 120, 180, 240, respectively. These values agreewithin uncertainty and show no clear trend. Similar be- FIG. 18. (Color online) Free energy density f as a function ofdyon density with S = 12, ν = for three different systemsizes. havior was seen for the SU (2) results [17], which studiedfinite volume effects in more detail and showed good re-sults for a maximum ensemble size of N D = 88. Weconclude that the number of dyons N D = 120 used inthis work is sufficiently large to determine the thermody-namic properties of the dyon ensemble. Appendix C: Topological charge fluctuations Here we present some general discussion of topologicalsusceptibility measurements in our setting and on thelattice.Suppose a lattice of 4-volume V is used, and the in-stantons with average density d I populate it randomly .This leads to Poisson distribution with mean number (cid:104) n I (cid:105) = d I V . The same is assumed for antiinstantons¯ I , and, as is well known, the susceptibility is simply thetotal density χ = (cid:104) ( n I − n ¯ I ) (cid:105) V = d I + d ¯ I (C1)In instanton liquid simulations the numbers N I , N ¯ I arefixed, and fluctuations appear only if one uses the sub-box as proposed in Ref. [49] for studies of correlationsin the instanton ensemble. Let the sub-box have vol-ume fraction f ≡ v /V . So, if dyons are also placed randomly , their distribution is binomial B ( f, N, n ) = N ! n !( N − n )! f n (1 − f ) N − n (C2)7In the limit when all volumes – and thus numbers in-volved – are large, standard statistical mechanics text-book analysis leads to the conclusion that both Poissonand binomial distributions lead to Gaussian fluctuationsaround their corresponding mean. Standard use of Stir-ling formula log ( n !) = nlog ( n/e ) and n = f N + δ leadsto B ( δ ) ∼ exp (cid:2) − δ f (1 − f ) N (cid:3) (C3)or (cid:104) δ (cid:105) = f (1 − f ) N . The topological susceptibility ofsub-box is then χ ( f ) = f (1 − f )( N I + N ¯ I ) f V = (1 − f ) χ P (C4)where χ P is that for unrestricted Poisson distributiondefined above. Indeed, the fluctuation must vanish at f → 1, but for small subsystem f → absolute value of the topo-logical susceptibility we measure to those on the lattice,one need to multiply the former by 1 / (1 − f ) = 2.Let us now calculate χ for randomly placed instanton-dyons. Their total numbers in our setting are also fixed, and we always keep N L = N ¯ L , N M = N ¯ M , N M = N ¯ M which ensures that the total magnetic and topologicalcharges of the whole system to be zero. The topologicalcharge in sub-box is Q = (1 − ν )( n L − n ¯ L ) + ν ( n M − n ¯ M ) + ν ( n M − n ¯ M )(C5)It is integer in two limiting cases, confining holonomy ν = 1 / T < T c , and trivial holonomy, ν = 0 at large T . The unrestricted topological susceptibility is then χ = (cid:0) (1 − ν ) ( N L + N ¯ L ) + ν ( N M + N M + N ¯ M + N ¯ M ) (cid:1) V (C6)In the confining case, with ν = 1 / N d , it is 6 N d / V . Comparing it to the casein which any triplet L, M , M is tightly coupled into in-stantons with the same density N I = N d , one finds that χ ( random dyons ) = 13 χ ( f ully correlated dyons ) (C7)Indeed, the clustering of dyons into instantons increasesthe fluctuations. Conversely, the clustering of dyons into D ¯ D pairs (as indeed happens at higher temperatures)decreases the fluctuations. [1] S. Mandelstam, Phys. Rept. , 245 (1976).[2] G. Parisi, Phys. Rev. D , 970 (1975).[3] K.-I. Kondo, S. Kato, A. Shibata, and T. Shinohara, Phys. Rept. , 1 (2015), arXiv:1409.1599 [hep-th].[4] G. S. Bali, in (1998) pp.17–36, arXiv:hep-ph/9809351.[5] A. D’Alessandro, M. D’Elia, and E. V. Shuryak, Phys. Rev. D , 094501 (2010), arXiv:1002.4161 [hep-lat].[6] J. Liao and E. Shuryak, Phys. Rev. C , 054907 (2007), arXiv:hep-ph/0611131.[7] A. A. Belavin, A. M. Polyakov, A. S. Schwartz, and Y. S. Tyupkin, Phys. Lett. B , 85 (1975).[8] E. V. Shuryak, Nucl. Phys. B , 93 (1982).[9] T. Sch¨afer and E. V. Shuryak, Rev. Mod. Phys. , 323 (1998), arXiv:hep-ph/9610451.[10] T. C. Kraan and P. van Baal, Phys. Lett. B , 389 (1998), arXiv:hep-th/9806034.[11] K.-M. Lee and C.-h. Lu, Phys. Rev. D , 025011 (1998), arXiv:hep-th/9802108.[12] N. Dorey, JHEP , 008 (2001), arXiv:hep-th/0010115.[13] A. Ramamurti, E. Shuryak, and I. Zahed, Phys. Rev. D , 114028 (2018), arXiv:1802.10509 [hep-ph].[14] E. Shuryak and T. Sulejmanpasic, Phys. Lett. B , 257 (2013), arXiv:1305.0796 [hep-ph].[15] Y. Liu, E. Shuryak, and I. Zahed, Phys. Rev. D , 085006 (2015), arXiv:1503.03058 [hep-ph].[16] R. Larsen and E. Shuryak, Phys. Rev. D , 094022 (2015), arXiv:1504.03341 [hep-ph].[17] M. A. Lopez-Ruiz, Y. Jiang, and J. Liao, Phys. Rev. D , 054026 (2018), arXiv:1611.02539 [hep-ph].[18] Y. Liu, E. Shuryak, and I. Zahed, Phys. Rev. D , 085007 (2015), arXiv:1503.09148 [hep-ph].[19] R. Larsen and E. Shuryak, Phys. Rev. D , 054029 (2016), arXiv:1511.02237 [hep-ph].[20] R. N. Larsen, S. Sharma, and E. Shuryak, Phys. Lett. B , 14 (2019), arXiv:1811.07914 [hep-lat].[21] O. Kaczmarek, F. Karsch, P. Petreczky, and F. Zantow, Phys. Lett. B , 41 (2002), arXiv:hep-lat/0207002.[22] D. J. Gross, R. D. Pisarski, and L. G. Yaffe, Rev. Mod. Phys. , 43 (1981).[23] D. Diakonov, N. Gromov, V. Petrov, and S. Slizovskiy, Phys. Rev. D , 036003 (2004), arXiv:hep-th/0404042.[24] M. Unsal and L. G. Yaffe, Phys. Rev. D , 065035 (2008), arXiv:0803.0344 [hep-th].[25] R. Larsen and E. Shuryak, Nucl. Phys. A , 110 (2016), arXiv:1408.6563 [hep-ph].[26] D. Diakonov, Nucl. Phys. B Proc. Suppl. , 5 (2009), arXiv:0906.2456 [hep-ph].[27] F. Bruckmann, S. Dinter, E.-M. Ilgenfritz, M. Muller-Preussker, and M. Wagner, Phys. Rev. D , 116007 (2009),arXiv:0903.3075 [hep-ph]. [28] V. G. Bornyakov, E. M. Ilgenfritz, and B. V. Martemyanov, Phys. Rev. D , 114510 (2020), arXiv:1908.08709 [hep-lat].[29] B. Alles, M. D’Elia, and A. Di Giacomo, Nucl. Phys. B , 281 (1997), [Erratum: Nucl.Phys.B 679, 397–399 (2004)],arXiv:hep-lat/9605013.[30] M. C`e, C. Consonni, G. P. Engel, and L. Giusti, Phys. Rev. D , 074502 (2015), arXiv:1506.06052 [hep-lat].[31] J. Frison, R. Kitano, H. Matsufuru, S. Mori, and N. Yamada, JHEP , 021 (2016), arXiv:1606.07175 [hep-lat].[32] G.-Y. Xiong, J.-B. Zhang, Y. Chen, C. Liu, Y.-B. Liu, and J.-P. Ma, Phys. Lett. B , 34 (2016), arXiv:1508.07704[hep-lat].[33] J. M. Pendlebury et al. , Phys. Rev. D , 092003 (2015), arXiv:1509.04411 [hep-ex].[34] F. K. Guo, R. Horsley, U. G. Meissner, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, A. Schiller, and J. M.Zanotti, Phys. Rev. Lett. , 062001 (2015), arXiv:1502.02295 [hep-lat].[35] E. Vicari and H. Panagopoulos, Phys. Rept. , 93 (2009), arXiv:0803.1593 [hep-th].[36] H. Panagopoulos and E. Vicari, JHEP , 119 (2011), arXiv:1109.6815 [hep-lat].[37] L. Giusti, S. Petrarca, and B. Taglienti, Phys. Rev. D , 094510 (2007), arXiv:0705.2352 [hep-th].[38] L. Del Debbio, H. Panagopoulos, and E. Vicari, JHEP , 044 (2002), arXiv:hep-th/0204125.[39] C. Bonati, M. D’Elia, H. Panagopoulos, and E. Vicari, PoS LATTICE2013 , 136 (2014), arXiv:1309.6059 [hep-lat].[40] C. Bonati, M. D’Elia, and A. Scapellato, Phys. Rev. D , 025028 (2016), arXiv:1512.01544 [hep-lat].[41] J. C. Myers and M. C. Ogilvie, Phys. Rev. D , 125030 (2008), arXiv:0707.1869 [hep-lat].[42] C. Bonati, M. Cardinali, and M. D’Elia, Phys. Rev. D , 054508 (2018), arXiv:1807.06558 [hep-lat].[43] C. Bonati, M. Cardinali, M. D’Elia, and F. Mazziotti, PoS LATTICE2019 , 084 (2019), arXiv:1912.12028 [hep-lat].[44] A. Athenodorou, M. Cardinali, and M. D’Elia, (2020), arXiv:2010.03618 [hep-lat].[45] E.-M. Ilgenfritz and E. V. Shuryak, Phys. Lett. B , 263 (1994), arXiv:hep-ph/9401285.[46] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lutgemeier, and B. Petersson, Nucl. Phys. B , 419(1996), arXiv:hep-lat/9602007.[47] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, JHEP , 056 (2012), arXiv:1204.6184 [hep-lat].[48] R. N. Larsen, S. Sharma, and E. Shuryak, Phys. Rev. D , 034501 (2020), arXiv:1912.09141 [hep-lat].[49] E. V. Shuryak and J. J. M. Verbaarschot, Phys. Rev. D52