Deconfinement Phase Transition in the $SU(3)$ Instanton-dyon Ensemble

Confinement remains one the most interesting and challenging nonperturbative phenomenon in non-Abelian gauge theories. Recent semiclassical (for SU(2)) and lattice (for QCD) studies have suggested that confinement arises from interactions of statistical ensembles of instanton-dyons with the 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 deconfinement phase transition. The properties of the ensemble, including correlations and topological susceptibility, are studied over a range of temperatures above and below $T_c$. Additionally, the dyon ensemble is studied in the Yang-Mills theory containing an extra trace-deformation term. It is shown that such a term can cause the theory to remain confined and even retain the same topological observables at high temperatures.


I. INTRODUCTION
Quantum Chromodynamics (QCD) is the quantum field theory describing the fundamental particles and forces that make up nuclear physics. While QCD is remarkably successful in describing nuclear physics, many phenomena remain beyond the scope of what can be studied analytically. Notably, nonperturbative phenomena such as confinement -the disappearance of quarks and gluons from the physical spectrum -is not completely understood. Confinement occurs not just in QCD, but in various Yang-Mills theories with or without quarks, making it clear that it emerges from the non-perturbative behavior of the gluons, rather than the quarks. Above certain critical temperature T c deconfinement takes place, and the QCD-like theories turn into a new form of matter, the Quark-Gluon Plasma (QGP).
Historically the first mechanism of the deconfinement transition was a 'dual superconductor' model [1][2][3]. At T < T c the chromoelectrically-charged quarks and gluons are connected by QCD flux tubes, dual to magnetic flux tubes in superconductor. With the advent of lattice gauge theories many aspects of this scenario were put to the test. In particular, the profile of the QCD flux tubes [4] was found to agree well with dual superconductor model. Monopoles were observed and found to rotate around these flux tubes, as expected. Bose-Einstein condensation of monopoles was detected and its critical temperature was shown to coincide with T c [5]. A high density of monopoles was found to be responsible for unusual kinetic properties of QGP [6].
Euclidean formulation of the gauge theory lead to discovery of 4D topological solitons known as BPST instantons [7]. A model of their ensemble, the Instanton Liquid Model (ILM) [8], has explained how instantons generate chiral symmetry breaking. As certain extrema of the path integral over gauge configurations, they form a basis for semiclassical theory, consistently including fluctuations around classical fields. Furthermore, one can study the interaction between instantons in their statistical ensemble: those studies explained behavior of correlation functions of various mesonic and baryonic currents, for review see e.g. Ref. [9]. Yet the instanton theory has not reproduced confinement.
Euclidean formulation of finite temperature QCD naturally led to a nonzero value of the Polyakov loop P = 0 as a signature of deconfinement. Note that throughout this paper we use P as shorthand for 1 3 T r[P ( x)] . It can be interpreted as the nonzero vacuum expectation value (VEV) of the time component of the gauge field A 0 , also known as nonzero holonomy. The natural question was then how to deform the instanton configurations in a way consistent with nonzero holonomy in the bulk. It was answered in Refs. [10,11], who discovered that instantons dissolve into N c (number of colors) constituent solitons, called instanton-dyons (or instanton-monopoles). Like original instantons, they are (anti)selfdual, so their actions and topological charges are equal. But, unlike instantons, their actions and topological charges are not quantized to integers; standard index theorems are avoided because instanton-dyons have magnetic charges and therefore are still connected by Dirac strings.
It was then realized that instanton-dyons provide a very valuable bridge between the theory of monopoles and instantons, providing a way to explain both confinement and chiral symmetry breaking in a single setting. Unlike monopoles, the instanton-dyons are semiclassical objects, allowing for the construction of a consistent theory of an interacting ensemble. Last but not least, the statistical sum in terms of instanton-dyons is "Poisson dual" to that based on monopoles, see Refs. [12,13].
Semiclassical approaches to finite-T gauge theories, with and without quarks, have lately been subject of multiple studies. Confinement in this theory is due to the back reaction of the dyon ensemble on the Polyakov loop, 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 semiclassical objects, unlike the QCD monopoles [13], with the actions S dyons ∼ 1/g 2 ∼ log(T /Λ) growing with temperature. Therefore their densities are suppressed at high T arXiv:2102.11321v1 [hep-ph] 22 Feb 2021 as an inverse power of T . At temperatures comparable to the critical one, T ∼ T c , numerically S dyons / ≈ 4. The inverse of it is the small parameter of our semiclassical expansion.
Although in this work we study pure gauge theory, rather than QCD-like theories with light quarks, let us mention that the instanton-dyon ensemble also describes the breaking of chiral symmetry. Numerical studies of chiral phase transition can be found in Refs. [18,19]. Using light quark Dirac eigenstates, one can identify the individual dyons inside lattice configurations from largescale QCD simulations. As shown recently [20], the lowest Dirac states are indeed generated by the instantondyons, in the form remarkably insensitive to large density of perturbative quarks and gluons in which they are immersed.
In this work we extend the studies [16,17] to the pure SU (3) gauge theory. Instead of two species of instantondyons, we now deal with three. It is well-known from lattice studies [21] that the pure SU (3) gauge theory (and all pure SU (N c ) theories with N c ≥ 3) possesses a first-order phase transition, rather than the secondorder transition seen in SU (2). So the question we address below is whether this, and other, features of pure SU (3) gauge theory can or cannot be reproduced by the instanton-dyon ensembles.
This paper is structured as follows: The dyon interactions and partition function are discussed in Section II. Section III lays out some technical details of the simulation as well as the data analysis performed. The physical results, including correlation functions and the temperature dependence of parameters such as the Polyakov loop are shown in Section IV. Finally, in Section V, the deconfinement transition of the dyon ensemble is studied in the trace-deformed Yang-Mills theory.

A. Dyons and holonomy
The Polyakov loop is an order parameter of the deconfinement phase transition where T a is the color generator in the fundamental representation. Its VEV, P is a unitary matrix with phases for eigenvalues, so A 4 can be defined as A 4 = 2πT diag(µ 1 , µ 2 , µ 3 ). These phases are related to the holonomies of the dyons by ν i = µ i+1 − µ i , where µ i+1 > µ i and µ 4 = µ 1 . The connection between the Polyakov loop and confinement can be seen through its relation to the free energy of a static quark Clearly, P = 0 is the confining holonomy, corresponding to F q → ∞, removing massive quarks from the spectrum completely. By combining the previous definitions of the holonomies and phases, the relationship between the holonomy and the average Polyakov line is seen to be From this it is clear that ν = 1/3 is confining ( P = 0),while ν → 0 produces the "trivial" Polyakov loop ( P = 1).
For SU (N c ) gauge theories, there are N c types of dyons: N c − 1 of them are called M i -type dyons, they correspond to the maximal diagonal subgroup. One more type is called the "twisted" (that is, 4-time-dependent) or L-type dyon. There are corresponding antidyons as well giving in total 2N c species of dyons. Note that we will use species to refer to one of the six dyons and type to refer to one of the three pairs of a dyon and its antidyon (i.e. an M 1 -type dyon refers to both M 1 andM 1 dyon species). The action of each individual dyon of type i is denoted by S 0 ν i , where ν i are the so-called holonomy parameters. Those are also called "fractions of the holonomy circle" because they satisfy the sum rule Nc i ν i = 1 and generically have N c − 1 independent parameters. These holonomy parameters determine not only the actions but also the core sizes of the dyons r i ∼ 1/(2πν i T ). In SU (N c ) theories, (anti)instantons are comprised of one of each of the types of (anti)dyons, although this is not necessarily the cause for other gauge groups. For example, in SU (3), the instanton I = M 1 LM 2 and the antiinstantonĪ =M 1LM2 .
The main parameter defining the properties of the dyons is the classical instanton action S 0 = 8π 2 g 2 ∼ log(T ) containing the coupling at the temperature-dependent scale. While this parameter is independent of the holonomy, the actions of individual dyons of type i are S i = S 0 ν i . In the SU (2) case, one defines one holonomy parameter ν ∈ [0, 1], the M dyon having action S M = S 0 ν and the L dyon having the conjugate action S L = S 0 (1 − ν). In the SU (3) case there are two diagonal Gell-Mann matrices, λ 3 and λ 8 and in principle two holonomy parameters. However, the average Polyakov loop P is a gauge-invariant, and thus physical, object. Restricting it to be real enforces ν M 1 = ν M 2 reducing the holonomy to just a single parameter for SU (3) with ν ∈ [0, 1/2]. The dyon actions expressed in it are This symmetry between the M 1 -and M 2 -type dyons means that they have equal densities as well, greatly reducing the space of parameters that needs to be considered. e interaction between thermal gluons and the holonomy. This generates the Gross-Pisarski-Ya↵e which appears in the an exponent in the partition function. For SU (3), this is given by shown with factors divided out so that it has units of free energy density. This potential disfavors minimum at trivial holonomy ⌫ = 0, and a maximum at the confining holonomy ⌫ = 1 3 . tion to the partition function consists of two parts, the contributions of the dyons in the absence ich can be expressed analytically directly from the input parameters, and Z np the contributions on, whose calculation is the subject of the simulation of the partition function, performed by ms in this work. The total partition function is their product Z = Z GP Y Z 0 Z np . ht for a single instanton was explicitly calculated in Ref. [23] for SU (2). Taking the dilute limit, ction between the dyons, this gives this factors into the weight for each individual dyon: Each dyon species contributes a factor of two holonomy terms stem from the holonomies of the individual M-and L-type dyons. A factor t for each dyon type. This is to remove a constant term that appears remains when taking the nowing this, it is easy to construct the instanton weight in SU (3) and extend it to an arbitrary summation. The partition function, assuming an equal number of dyons and antidyons, is t of an individual dyon with holonomy ⌫, the limit V ! 1 and assume equal densities of M 1 -and M 2 -type dyons. The free energy

B. The partition function and dyon interactions
A complete calculation of the dyons' free energy requires the construction of the dyonic partition function. We start first with effects that are not induced by the dyons' non-perturbative interactions. In the absence of all dyonic effects, there is a perturbative interaction between thermal gluons and the holonomy. This generates the Gross-Pisarski-Yaffe potential V GP Y [22], which appears in the an exponent in the partition function. For SU (3), this is given by Here the potential is shown with factors divided out so that it has units of free energy density. This potential disfavors confinement, having a minimum at the trivial holonomy ν = 0, and a maximum at the confining holonomy ν = 1 3 . The dyon contribution to the partition function consists of two parts, the contributions of the dyons in the absence of interactions Z 0 , which can be expressed analytically directly from the input parameters, and Z np the contributions of the dyon interaction, whose calculation is the subject of the simulation of the partition function, performed by Monte-Carlo algorithms in this work. The total partition function is their product Z = Z 0 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 It is easy to see how this factors into the weight for each individual dyon: Each dyon species contributes a factor of √ ΛS 0 e − √ Soν and the two holonomy terms stem from the holonomies of the individual M-and L-type dyons. A factor of 4πν i is divided out for each dyon type. This is to remove a constant term that appears remains when taking the dilute limit in Z np . Knowing this, it is easy to construct the instanton weight in SU (3) and extend it to an arbitrary number of dyons by summation. The partition function, assuming an equal number of dyons and antidyons, is where d ν is the weight of an individual dyon with holonomy ν, Now we may take the limit V → ∞ and assume equal densities of M 1 -and M 2 -type dyons. Additionally we resolve the factorial terms with Stirling's approximation carried out to three terms ln N ! ≈ N ln N − N + ln( √ 2πN ). The free energy F = − ln(Z) is given by the following expression where ∆f is the free energy density stemming from the interactions of the dyons. If the dyons have classical binary interactions ∆S class and a volume metric G, their contributions to the partition function and the free energy density are The set of parameters that minimizes the free energy density corresponds to the physical dyon ensemble in the infinite volume limit. This elucidates the main procedure of this work: to first compute the free energy density of the ensemble for a wide range of parameters, and then locate their values which minimize it. The classical interactions between dyons and antidyons are, at distances exceeding the dyon cores, asymptotically Coulomb-like. For generic SU (N c ) theories, the interactions are given in Ref. [24]. At shorter length scales, the interaction is modified. The interaction between dyons and antidyons of the same type was studied in detail and parameterized by Larsen and Shuryak [25]. In the previous SU (2) model, this interaction was used for dyon-antidyon pairs of the same type, while dyons and antidyons of different types had purely Coulomb-like interactions. In this work we use a single modified parameterization for all dyon-antidyon pairs. The parameterization used is given by where C dd is a coefficient with value 2 for pairs of the same type and −1 for pairs of different types. Compared with the parametrization used in SU (2), two modifications have been made: the substitution ν i → √ ν i ν j was made to accommodate pairs of dyons with different holonomies (e.g. M 1L ), and the coefficient in front of the exponential term has been reduced from the original values of 3.264. These interactions are used for distances greater than the core size r 0 = x 0 /(2πνT ). Dyons pairs of the same type (regardless of duality) experience a repulsive core. While the core potential has not been studied in detail it reasonably described by The two parameters in this interaction, x 0 and V 0 are chosen on a phenomenological basis and should be subject to constraints from appropriate lattice data when possible. See Appendix B for a discussion of their values and the effects of changing them. Additionally, the dyons experience an effective potential from the fluctuation determinant of the instanton. This effect leads to the Diakonov determinant [26] of the metric of the space of dyons' collective variables where r i,m is the position of the i'th dyon of type m. This metric only accounts for the selfdual dyons. An equivalent metricḠ is used for the antidyons as well.
The metric is true for dyons of different species at any distance, but only true at larger distances for dyons of the same species. We modify the terms with a cutoff r → r 2 + (3/2πT ) 2 so that the diagonal entries go to 0 rather than −∞ when ν = 1 3 . The effective potential from this metric comes as The diagonal of this metric does not vanish in the dilute limit, instead going to i 4πν i . This has been accounted for by modifying the individual dyon weight in Eq. (8).
It is observed that in the case the densities of different dyon types are unequal, then in the infinite volume limit, the sum diverges. We therefore regulate all Coulomb terms in both the classical and one-loop potentials with the dimensionless Debye mass r → re M D rT . Like the core parameters, the value of the Debye mass is a choice that should be subject to improvement. Now let us give a brief qualitative discussion of how such interactions should generate confinement. The main interaction driving confinement is the repulsive cores of the dyons. Because the volume of the cores goes as 1/ν 3 i , this interaction disfavors any type of dyon becoming large and thus favors the confined phase ν = 1 3 where all dyons are the same size. The long-distance classical Coulomb interactions disfavor confinement however. Because dyons repel anti-dyons of different types but attract antidyons of the same type, these interactions favor ensembles with many dyons of the same kind -a large value of n M /n L -rather than equal numbers of all species. At larger values of n M /n L the entropy portion of the partition function Z 0 more strongly favors small values of ν, driving the system to the deconfined phase. The one-loop Coulomb-like interactions have the opposite sign and thus opposite effect. These one-loop interactions are suppressed by a factor of S 0 and are only comparable to the classical interactions at low T .
The nonperturbative interactions grow weak compared to the perturbative portions of the partition function as T increases. At low temperature we have an ensemble where the cores and one-loop interactions dominate, binding the dyons into instantons and at high temperatures we have a gas of mostly M iMi pairs bound by their classical attraction and the entropy and Gross-Pisarski-Yaffe terms have shifted ν to a lower value.

III. THE SIMULATION SETTING AND DATA ANALYSIS
Numerical integration over the dyons' coordinates is performed by Monte-Carlo techniques. The position of each dyon is updated and then accepted or rejected according to the standard Metropolis algorithm. Once this is done for each dyon five times (five updates is approximately the autocorrelation time of the ensemble), the configuration is sampled and its properties are computed. The free energy is obtained via the usual integration over auxiliary parameter λ: where this nonperturbative free energy is added to the perturbative contribution from Eq. 9. Integration over the dummy parameter is performed in 10 equal steps λ = 0.1, ... 1, with the average energy at each step being computed from 2000 configurations. Because the contribution at small λ is large, the value of the average interaction at λ = 0 is obtained from a linear fit of the nearest points and included in the numerical integration. With all of this, the statistical uncertainty in the computations of f are generally at the percent level. The goal of these simulations is to calculate the free energy density f of ensembles with three main input parameters, the holonomy ν and the dyon densities n M , n L over a range of temperatures near T c . All simulations are performed with a fixed number of dyons N D = 120 or 122 (60 or 61 dyons and 60 or 61 antidyons). For the six species of dyons in SU (3), N D = 4N M + 2N L . The system is studied on flat geometry with finite density imposed by a periodic box surrounded by 26 image boxes, as was done for SU (2) in Ref. [17]. The densities of the dyons, n M and n L are controlled by varying the size of the 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 are used to allow for more values of N M /N L to be used.
There are two technical details of the simulation that must be mentioned. The first is the Diakonov determinant det G, which is not positive definite for randomlyplaced dyons. It has been shown by Bruckmann et. al. [27] that in fact, at reasonable densities, nearly all random configurations have at least one negative eigenvalue λ i . This is handled by a modification of the effective potential where V max is a large value (we choose V max = 100) which serves to reject all Metropolis updates into the region of negative eigenvalues. The choice of V max is arbitrary so long at it is large enough to ensure that no configuration with negative eigenvalues is accepted during the course of the simulation. The second is an approximation which must be made due to computational resources. The total system, the main box and all its nearest images, contains 27N D dyons. Its 2D version is shown in Fig. 2. The most demanding step in the simulation is the computation of det G which is an O(N 3 ) process. Performing this computation at every update step of the Metropolis algorithm is numerically costly, given the large number of parameters for which the free energy must be calculated. Instead, the change in action associated with each update is computed considering only the dyons in a box of the same size as the simulation box centered at the position of the dyon being updated. This approximation 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 depen- in the simulation is the computation of det G which is an O(N 3 ) process. Performing update step of the Metropolis algorithm is numerically costly, given the large number of pa energy must be calculated. Instead, the change in action associated with each update is the dyons in a box of the same size as the simulation box centered at the position of the approximation is well-justified by the fact that all long-range interactions are regulated by Once the free energy density of each set of parameters is determined, a fit is performed to of the free energy density near the minimum. For each value of the action S 0 there i input parameters with coordinates (⌫, n M , n L ). The points near a minimum are fit to a coordinates. From the fit, the minimum value of f and its location in the coordinate spac value of S 0 tested, all local minima are identified and the points near them are fit as describ minimum value of f is the global minimum. For cases where a minimum is located in the ⌫ = 1 3 and n M = n L , these equalities are assumed exact and the minimum is fit to a quad 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 (ii) determine the global minimum for each temperature via fitting to determine the physica ⌫(T ), n M (T ), n L (T ), (iii) study the ensembles with the physical values with increased statistics to compute c susceptibility. (iv) With these increased-statistics runs, we generate 50,000[?] configurations with the ph them ( = 1, we no longer need to integrate over it) using the fitted values of the inputs for In order to more accurately reproduce the fitted values of both dyon densities simultaneou allowed to vary more with 126  N D  164. dence of the free energy density near the minimum. For each value of the action S 0 there is a 3-dimensional space of input parameters with coordinates (ν, n M , n L ). The points near a minimum are fit to a quadratic function in these coordinates. From the fit, the minimum value of f and its location in the coordinate space is determined. For every value of S 0 tested, all local minima are identified and the points near them are fit as described. The fit with the lowest minimum value of f is the global minimum. For cases where a minimum is located in the in the confined phase with ν = 1 3 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 topolog-ical charge distributions. (iv) With these increased-statistics runs, we generate 80,000 configurations with the physical interactions between them (λ = 1, we no longer need to integrate over it) using the fitted values of the inputs for each temperature studied.

IV. RESULTS
A. Structure of the holonomy potential and the phase transition It is expected that the pure SU (3) gauge theory possesses a first-order phase transition at some critical temperature T c . The dyons interactions generate the holonomy potential which has a global minimum in the confined phase at T < T c and then a jump to a new global minimum in the deconfined phase at T > T c . The first goal of this work was to check for the existence of such a phase transition. Parameters of the interaction were chosen so this occurs around S 0 ∼ 12 (See Appendix B for more details). For each value of the temperature, which is related to S 0 by Eq. (A1), the shape of the holonomy potential was computed and the minimum was determined. The value of the holonomy and dyon densities at which this minimum occurs are the physical values the ensemble takes on at that temperature.
The instanton action S 0 was varied in steps of 1 from S 0 = 8 -21. It was observed that the phase transition occurs between 13 < S 0 < 14 and additional simulations were performed at S 0 = 12.75, 13.25, 13.5, and 13.75 to locate the critical temperature more accurately. A linear fit was performed to the free energy density with the nearest points on both sides of the transition to determine that the phase boundary is located at S 0 = 13.18. With this, the action can now be rewritten in terms of T /T c . The range of temperatures studied here is then 0.62T c < T < 2.04T c . In physical units, this is approximately 160 MeV < T < 530 MeV.
This model of the dyon interactions relies on a semiclassical expansion. An important consideration is then, whether or not the main semiclassical quantity, the action of the dyons S i = S 0 ν i , remains sufficiently large over this temperature range. In the confined phase, all dyons have action S 0 /3, meaning there is clearly a lower bound where the action becomes small and we choose S 0 = 8, S i = 2.6 as a reasonable lower bound. In the deconfined phase, as the instanton action grows, the holonomy decreases making the action of the M i -type dyons the limiting factor. The action of these dyons remains nearly constant at S M 3 up to 2T c . Thus, the dyons remain sufficiently semiclassical throughout the entire range of temperatures studied.
By performing similar linear fits, all parameters of the ensemble can be determined on both sides of the phase transition.  At temperatures below T c the ensemble is in the confined phase with ν = 1/3 and n M = n L , as can be seen in Fig. 3 (left). At densities below the physical one, the minimum shifts to the left as the nonperturbative interactions become weak compared to the perturbative contribution to the free energy. At higher densities the ensemble prefers to remain in the confined phase, but with a larger free energy minimum and curvature of the potential.
The deconfined phase has a similar structure, but with the global minimum occurring at ν < 1/3 and n M > n L . However at densities above that of the global minimum, the minima continue to move towards larger ν. At these densities, the repulsive core dominates and it becomes energetically favorable to make the many M i dyons smaller at the cost of making the few L dyons larger. It is possible that at densities higher than what were studied here, the ensemble may have a minimum at ν > 1/3.
These plots only show a slice of the full space of parameters explored for each value of S 0 for specific values of n M /n L . The structure of the first-order phase transition can be seen more clearly in Fig. 4. By considering the minimum free energy density selected from all combinations of n M and n L as a function of the holonomy, the two local minima -one in the confined phase and the other in the deconfined -are clearly visible. At the value of S 0 nearest to the critical value, the free energy at the two minima are nearly degenerate and the global minimum switches between the two as the temperature changes.
This structure is different from the that of SU (2). In SU (2), where the phase transition is second order, rather than having two degenerate minima, the holonomy potential flattens near T c (see e.g. Fig. 5 of Ref. [16]). This allows the minimum to quickly, but smoothly, shift from the confining holonomy to smaller values. Additionally, there is a ν ↔ 1 − ν symmetry not present in SU (3).

B. Temperature dependence of the parameters
The free energy density f , unlike other physical quantities, remains continuous across even a first-order phase transition, as we see in Fig. 5. Its derivative, however, may not. The free energy varies with temperature much more rapidly in the confined phase than the deconfined. The most important feature of the dyon ensemble for describing the deconfinement transition is the average Polyakov loop as a function of the temperature P (T ) . Below T c , the holonomy takes the confining value ν = 1/3, P = 0. At T c the value jumps to ∼ 0.4 and then continues to increase as T increases. The value of the average Polyakov loop above the phase transition shows qualitative agreement with the lattice data [21], but does not increase with temperature as quickly. A change to the parameters of the dyon interactions could improve the agreement with the lattice data.
Another set of properties of the ensemble are the densities of each dyon type, shown in Fig. 7. At T < T c , all densities are equal reflecting the fact that, in the confined phase, all dyons have equal statistical weights, core sizes, and masses and have symmetric interactions between them. At T > T c , the holonomy decreases discontinuously, causing a similar change in both dyon densities. In the case of the M -dyon density n M , the decrease in ν increases the statistical weight of the M i -type dyons, while simultaneously increasing the size of their repulsive cores. These competing effects result in what is seen in the ensemble: a small, O(5%), decrease in n M at the phase transition. The L-dyon density sees a more significant decrease due to the large decrease in statistical weight of the L dyons from the increase in 1 − 2ν.
Overall, the dyon densities decrease with temperature, consistent with the expected behavior of instantons. It is the dyon densities that set the upper limit on the temperature that can be studied in our ensemble. As temperature increases, the ratio n M /n L becomes larger and larger, and, at the largest values, the number of L dyons in 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 0 = 21) is a technical constraint, rather than a physical one; accurately probing higher S 0 would require a larger ensemble with better statistics.

C. Spatial correlations between dyons
It is useful to study the effects of the dyons' interactions by studying the spatial correlations between them. Such correlations will be useful for comparison to studies of the dyons in lattice configurations. The most straightforward characteristic of the spatial correlations are the correlation functions C ij (rT ) between two dyons of species i and j. In SU (3), because the instanton is comprised of three constituent dyons, it is useful to also define the 'instanton correlation function' C I (yT ), where y is the hyperdistance in the 6-dimensional space of Jacobi coordinates of the three dyons In both cases, the functions are normalized such that C ij , C I → 1 at large (hyper)distance and have had the angular factors divided out ((rT ) 2 for the two-dyon correlation functions and (yT ) 5 for the instanton correlation function). Many of the correlation functions are redundant so we need to only observe a subset of the correlation functions to understand the behavior of the ensemble.  In several channels we can see clearly the effects of the repulsive core at short ranges. In the DD channels, we can see positive correlation just beyond the cores. Changes in the cores can also be seen: in the M 1M1 channel, the core becomes larger and softer above T c , soft enough that many pairs can still be found at x < x 0 , while the cores of the L dyons becomes smaller and harder, reducing the correlation function to 0 at x < x 0 . The smaller core size also allows the attractive long-range interactions to become larger near to the core, enhancing the correlation found there.
Despite having no long-range interactions, there is small correlation seen in the same-species (M 1 M 1 and LL) channels beyond the core due to some mutual correlations to other dyons. Compared the the DD channels the anticorrelation at small r is even stronger due to additional repulsion from the Diakonov determinant. Some correlation functions for the dyons in pure SU (3) Yang-Mills theory at T /T c = 0.96 have been observed in lattice configurations [28]. We find qualitative agreement with these results in all channels except for dyons of the same species. Here some short-range correlation is observed among the lattice dyons, suggesting an absence of the repulsive core used in our model in this channel.

D. The vacuum angle and moments of the topological charge
Another set of important characteristics of the dyon ensemble are its topological properties. Non-abelian gauge configurations generically may possess some topo-logical charge Q with its density q(x) given by The (anti)instantons are topologically nontrivial objects with Q = ±1. The constituent dyons then each carry some fraction of this charge. Dyons of species i carry topological charge Q = ν i while antidyons have charge Q = −ν i . Because we only consider ensembles with equal numbers of dyons and antidyons, Q = 0 in any sub-volume of the box. The most prominent topological feature of the gauge theory, which has been extensively studied on the lattice [29][30][31], is the topological susceptibility χ = Q 2 /V 4 . Because the dyons in our semiclassical ensemble have well-defined positions and topological charges, measuring the average topological charge only requires summing over the charges of the dyons in a given sub-volume. This is done by splitting the main simulation box in half along each of the 3 axes and computing the total charge in the said half-boxes. This results in 3 independent measurements per configuration, meaning the average of any power of Q is computed from 240,000 measurements.
The temperature dependence of the topological susceptibility over the entire temperature range probed in this work is shown by blue points in Fig. 10. To make comparison to lattice data easier, we show the relative susceptibility, normalized all sets to their values just below the critical temperature. The issue of absolute comparison is further discussed in Appendix C.
With χ(T ), like P (T ) , we see a clear first-order phase transition in the data. Unlike with the Polyakov loop however, such an abrupt transition with a jump is not obvious from the lattice data points (although they of course do not contradict existence of a jump). The topological observables have been analyzed on the lattice in many works, we compare our dyon results to those in Refs. [29,32], which cover a similar temperature range as we have studied.  [29,32]. Lattice values 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 point size.
The (Euclidean) Lagrangians of generic SU (N ) gauge theories can be appended by an additional topological term where θ is the so-called vacuum angle and q(x) is the topological charge density as previously defined in Eq. (20). A non-zero θ explicitly breaks CP symmetry. In QCD it is known that |θ QCD | 10 −10 , as it is constrained by the upper bound on measurements of the neutron's electric dipole moment [33,34].
By expanding the free energy density f (θ) around θ = 0, its dependence on θ can be studied at small, but nonzero θ. It can be expanded as [35] where the coefficients b n are related to cumulants of the topological charge computed at θ = 0. We compute the first two terms, which are given explicitly by In the confined phase, we find that all b 2 values are consistent with b 2 being constant below T c . These seven values vary around 0.02 − 0.03 and have an average b 2 (T < T c ) = 0.026. Above T c , b 2 quickly drops and then remains approximately constant just above 0.01. This behavior is in disagreement with available lattice data. The T = 0 value of b 2 has been particularly well studied on the lattice [36][37][38] with values around b 2 (0) = −0.025. On the high-temperature end, lattice data [32,39] finds that b 2 approaches the value predicted by the Dilute Instanton Gas Approximation (DIGA), b 2 (T ) → −1/12.
Below T c , our values of b 2 are consistent with the magnitude of the T = 0 value predicted on the lattice, but with opposite sign. In the high temperature limit, we again see that our dyon model predicts positive values rather than negative, as well as being an order of magnitude smaller than the lattice results. Clearly, these small non-Gaussianities in the topological charge distribution are quite sensitive to the dyon interactions. Changes to the dyon interactions could improve the agreement with lattice data.
Finally, our values for b 4 are compatible with zero for all temperatures as is also observed on the lattice [36,40]. Upper limits from these lattice results constrain its value to be |b 4 (T = 0)| < 10 −3 . Trace-deformed gauge theories (TDGT) were introduced in Refs. [24,41], adding certain terms containing powers of the Polyakov line to the theory action, such as with a new parameter h. At T > T c , when in the undeformed theory the Polyakov loop is nonzero, the new term obtains an additional contribution and suppresses its value, shifting the theory back toward the confining holonomy. For large enough h the system returns to the P = 0 phase, which we will call the "reconfined phase". These theories have been studied in numerous lattice simulations, of which we will mention those of the Pisa group [42][43][44]. Their main findings are that for a specific temperature above T c at large enough h > h * ∼ 1/3, at which P returns to zero, there is no more dependence on h, and the reconfined phase is remarkably similar to the original (low temperature) confined phase. It was shown to possess the same spectrum, topological observables, and even spectrum of periodic strings (torelons).
From the perspective of the instanton-dyon theory that we are developing, it is clear that the reconfined phase, like the confined one, corresponds to symmetric value of the holonomy ν = 1/3 and the same densities of all species L, M 1 , M 2 of the dyons. Yet the temperature scale, and thus the periodicity interval of the Euclidean time τ , is different. Therefore, the absolute scale of the dyon action parameter S 0 = 8π 2 /g 2 (T ) must increase, making the system more dilute. And yet, the topological susceptibility χ is the same as at T = 0 despite a significantly reduced dyon density! It is possible to study the instanton-dyon ensemble in a theory with arbitrary h by the addition of a deformation term to the free energy density ∆f def = h P 2 . Because this new term depends only on the value of the holonomy, one does not need to perform new MC simulations to study the effect of changing h; one needs only to add the new term to the un-deformed results and perform new fits to find the minima.
The first questions to ask are how does changing h affect the value of the Polyakov loop and at what value of h does the system become reconfined? Fig. 12 shows the temperature-dependence of the Polyakov loop for a few different values of h. As expected, the larger the trace deformation, the more suppressed the values of P 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 confining value. The critical temperature of the trace-deformed theories are determined via the same method as the original theory: the intersection of linear fits to f on both sides pf the transition. At this point we are limited to determining T c for h ≤ 40, as above this value the the Polyakov loop remains very close to the confining value and the minimum value of the holonomy fit become closer to ν = 1 3 than the holonomy step size used for the simulations, meaning that ν = 1 3 is compatible with our uncertainties for all temperatures studied. In order to accurately determine the critical temperature for large h, smaller steps in both S 0 and ν are necessary. This also makes determining the nature of the phase transition inconclusive at larger h. At small h, it is clear to see that the phase transition remains first order, while at larger h (above ∼ 10) we do not have sufficient resolution in ν to determine whether the holonomy potential continues to have two distinct minima or becomes a smooth crossover. While we are currently limited in how high of a critical temperature we can probe (up to ∼ 1.3T c (0) as per Fig.  13), the theory may have non-trivial behavior at higher temperatures. In the limit that T → ∞ the density of dyons goes to zero and only two terms in the free energy density remain: the (deconfinement-favoring) Gross-Pisarski-Yaffe potential and the (confinement-favoring) trace-deformation term. Clearly which term dominates depends on the value of h. At h > 5π 2 /18 2.74 the confining holonomy is the global minimum.
This leads to the conclusion that there are three distinct regimes of the TDGT: i) h < 5π 2 /18: The theory is confined at low temperatures and deconfined at high temperatures. ii) 5π 2 /18 < h < h max : The theory is confined at low temperatures, deconfined in some intermediate region and 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 temperatures to directly see the second phase transition in regime (ii), the last value of P for h = 40 in Fig. 12 suggests that P may be starting to decrease and return to zero.

B. Topological observables of the reconfined theory
One of the most interesting features of the deformed gauge theory observed on the lattice [42][43][44] is that when the trace deformation is large enough to return the system to the confining holonomy, the topological observables, namely χ and b 2 , also return to the same values they 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 as was done in the lattice studies, we choose an arbitrarilylarge value of h and vary temperature. As mentioned in the previous section, for some large value of h, such that h > h max , the system should remain in the confined phase at all temperatures. In our dyon model, we can study this scenario, which we call the 'maximallydeformed' theory, by simply demanding that ν = 1 3 and determining the dyon densities that minimize the free energy density with this constraint. Because our individual simulations have fixed holonomy value ν = 1 3 , P = 0 the trace-deformation term does not contribute to the free energy density of these simulations and doesn't need to be considered in the fits. The question then is whether or not the dyon ensemble, like the lattice theory, possesses the same topological observables in the maximally-deformed theory as in the original confined phase. In the case of χ, our dyon model did not predict a constant value but an approximately linearly-decreasing one. We observe analogous behavior to what is seen on the lattice; in the maximally-deformed Topological susceptibility of the maximallydeformed theory as a function of temperature. Error bars are smaller than the point size.
theory the abrupt jump in χ disappears and it displays the same temperature dependence above T c (0) as below. Comparing the dyon density (Fig. 14) and the topological susceptibility (Fig. 15), one can see that the decrease in χ is mostly driven by the decreasing dyon density. Although n i decreases by a factor of ∼ 4 while χ decreases by a factor of ∼ 7 over the temperature range we study, meaning the value of Q 2 decreases by a factor of ∼ 2 as we indeed observe.
The behavior of b 2 is more inconclusive. It remains constant up to about 1.1T c (0) and then decreases quickly and is compatible with zero at all temperatures above that point. Again, it remains positive at all temperatures. It should be noted that the error bars on the b 2 measurements are still quite large and the trend could look quite different with improved statistics.
In order to determine what drives this behavior, it is useful to look at some correlation functions in this deformed theory. It is clear that the ensemble has some non-trivial temperature dependence even though it remains in the confined phase. The two main channels with positive correlation are, just like the un-deformed theory, the DD channels and the instanton (M 1 M 2 L) channel. As the temperature is increased, the ensemble more strongly prefers DD pairs rather than (anti)instanton as can be seen in Fig. 17. This is due to the increase in the strength of the DD classical attraction. Eventually, these pair correlations grow strong enough to completely kill correlation in the instanton channel.
In the original theory, the system prefers binding into instantons in the confined phase and then transitions into preferring M 1M1 -and M 2M2 pairs at high temperatures; the density of LL is becoming small at this point. However, 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 behavior not seen in the original theory. At high temperature the system becomes a three-component gas of DD pairs, with all types having equal densities and core sizes. Because the theory remains confined, the core sizes never become large and the correlation between dyons in DD pairs becomes very large. The topological observables are very sensitive to this binding , as the (anti)instantons and DD pairs have topological charges Q = ±1 and Q = 0, respectively. The shift towards the latter leads to a neutralization of the topological charge and is responsible for driving the values of χ and b 2 down at higher temperatures.
Looking to the future, one can see that this behavior has interesting implications when quarks are added to the theory. It is known that chiral symmetry restoration is directly related to the quark zero modes on the instantons. In particular, the quark interactions drive the system to form instanton-antiinstanton 'molecules'. This neutralizes the fluctuations of the topological charge and shifts the Dirac eigenvalues away from zero at high temperatures, restoring chiral symmetry [45]. In fact in this work we already observe DD pairing at higher temperatures. We see here that such correlations continue to grow, even when the theory remains confined due to the deformation term in the action. This suggests that in the trace-deformed theory with quarks one may see chiral symmetry restoration occur inside the confined phase.

VI. SUMMARY AND DISCUSSION
This paper reports the first direct numerical simulations of the ensemble of instanton-dyons in the pure SU (3) Yang-Mills theory. We use classical and semiclassical one-loop measures which have derived the interdyon interactions. We have performed Monte-Carlo simulations of an ensemble of dyons in a 3D cube with periodic boundary conditions over a range of temperatures, densities, and holonomy values.
Our main objective is to see whether the semiclassical ensemble of instanton-dyons does or does not reproduce the deconfinement transition, which is, in this case, of the first order. We indeed find that there are two distinct phases, the confined (with free energy minimal at holonomy value 1/3, Fig. 3 (left)) and the deconfined (with free energy minimal elsewhere, Fig. 3 (right)). We were able to map out the holonomy potential and see explicitly the two nearly-degenerate minima of the potential representing the transition between the two phases (red squares in Fig. 4) and study various properties near the critical temperature.
Comparison between our semiclassical model and lattice data for VEV of the Polyakov loop P (T ) is shown in Fig. 6. Remarkably, the jump of it at T c is reproduced rather precisely, with some moderate deviations at higher temperatures. As shown in Fig. 7, this large jump results in L-dyon suppression in the deconfined phase, while the density of M i -dyons has a rather smooth temperature dependence.
The pair-wise spatial correlations between various dyons are shown in Fig. 8. Overall, they are in qualitative agreement with expectations. Note also, that for a system with a strong first order transition, one does not see drastic changes in the correlations, except for overall change 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 typical distribution is shown in Fig. 9, at first sight having just a normal Gaussian form. Its width -the topological susceptibility χ -is compared to lattice data in Fig.  10. The magnitude of the jump and behavior in the deconfinement phase are found to be quite similar, but the dependence at T < T c is somewhat different. Our results for non-Gaussianity parameter b 2 are perhaps still too noisy. Let us note that these observables, as well as the dyon correlations, are very sensitive to the details of the dyon interactions, particularly the short-range classical interactions, which are the strongest. These interactions are at present the most poorly-understood aspect of the dyon model with some parts being analytically derived (e.g. the Diakonov determinant) and others being simply phenomenological (e.g. the repulsive core). In order to achieve a quantitative agreement between the dyons and lattice data, if possible, these interactions should be the subject of further study and improvement.
We also report the first study of trace-deformed gauge theory in the instanton-dyon setting. As expected, by increasing the coupling h of the deformation term, one indeed suppresses the jump and value of the Polyakov loop in the deconfined phase. The phase transition location is also pushed to higher values, see Fig. 13. Overall, like on the lattice, we see that increasing h leads to a return to the topological susceptibility as in the confining phase of the un-deformed theory. We however still see that this agreement is only partial, see for example the spatial correlations taken in the maximally-deformed confining theory.
Completing this summary of our results, we again emphasize that the semiclassical model used, which has many orders of magnitude fewer degrees of freedom than lattice gauge theory, is able to explain properties of the deconfinement phase transition of the SU (3) gauge theory quite well. Some quantities -like the jump in the Polyakov line -are reproduced precisely, others semiqualitatively, but we have not seen any serious disagreements. The double-trace deformation study also shows that both our model and lattice react to this addition to the action in a very similar way.
Thinking of the future, the next step will of course be the inclusion of dynamical quark flavors via their zero modes and interactions. Such work will allow not only for the study of confinement (which will now be expected to be a smooth crossover), but also to study the role of the dyons in chiral symmetry restoration.
were chosen. This increases the number of configurations to be tested, but allows for a more general, temperaturedependent Debye mass. Both methods have been shown to produce reasonable properties for the SU (2) dyon ensemble, so we have saved computational time by treating it 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.

Finite-size effects
In principle, equilibrium statistical mechanics is derived in the thermodynamic limit N D ,Ṽ 3 → ∞, which is not possible in numerical simulations. For such a finite numerical simulation, the main concern is whether the system size used is sufficiently large enough so that the measured observables can be reliably extrapolated to the thermodynamical limit. The simplest way to check these effects is to simulate larger systems with the same parameters and study the dependence of some observables on the size N D .
Here we discuss a representative example of the finitesize effects as seen in Fig. 18. In increasing the size of the system by up to a factor of 2, there is a noticeable change 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 stemming from Stirling's approximation, which has explicit V 3 -dependence. Increasing the system size decreases the contribution from the third term as it goes as 1/Ṽ 3 . Subtracting this term from the results reveals that the contribution from the dyon interactions increases with system size. Both of these effects are reduced in the deconfined phase where the density is lower.
The main factor in the finite-size effects regarding dyon interactions comes from dyons near the boundaries of the system. As long as a finite number of dyons and image boxes are used, there are dyons near the faces of the cube whose total short-range interactions are reduced. The fraction of dyons near the faces of the cube should scale as N −1/3 D and thus vanishes in the thermodynamic limit. Parameters that determine the effective range of interactions such as the core sizes and Debye mass also play an important role in determining how quickly the results 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 energy minimum decreases by ∼ 2% when doubling the number of dyons. A simple linear fit of f as a function of 1/N D gives an infinite-volume value of f = −4.519, a 4.5% decrease from the N D = 120 value. The input parameters of the ensemble -the holonomy and densities -depend only on the location of the minimum, which clearly varies much less. Fitting the curves to quadratic forms gives the densities n M = 0.607, 0.609, 0.607 for dyon numbers N D = 120, 180, 240, respectively. These values agree within uncertainty and show no clear trend. Similar be- havior was seen for the SU (2) results [17], which studied finite volume effects in more detail and showed good results for a maximum ensemble size of N D = 88. We conclude that the number of dyons N D = 120 used in this work is sufficiently large to determine the thermodynamic properties of the dyon ensemble.

Appendix C: Topological charge fluctuations
Here we present some general discussion of topological susceptibility measurements in our setting and on the lattice.
Suppose a lattice of 4-volume V 4 is used, and the instantons with average density d I populate it randomly. This leads to Poisson distribution with mean number n I = d I V 4 . The same is assumed for antiinstantons I, and, as is well known, the susceptibility is simply the total density In instanton liquid simulations the numbers N I , NĪ are fixed, and fluctuations appear only if one uses the subbox as proposed in Ref. [49] for studies of correlations in the instanton ensemble. Let the sub-box have volume fraction f ≡ v 4 /V 4 . So, if dyons are also placed randomly, their distribution is binomial In the limit when all volumes -and thus numbers involved -are large, standard statistical mechanics textbook analysis leads to the conclusion that both Poisson and binomial distributions lead to Gaussian fluctuations around their corresponding mean. Standard use of Stirling formula log(n!) = nlog(n/e) and n = f N + δ leads to or δ 2 = f (1 − f )N . The topological susceptibility of sub-box is then where χ P is that for unrestricted Poisson distribution defined above. Indeed, the fluctuation must vanish at f → 1, but for small subsystem f → 0 both become the same. So, in comparing the absolute value of the topological 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 instantondyons. Their total numbers in our setting are also fixed, and we always keep N L = NL, N M 1 = NM 1 , N M 2 = NM 2 which ensures that the total magnetic and topological charges of the whole system to be zero. The topological charge in sub-box is Q = (1 − 2ν)(n L − nL) + ν(n M 1 − nM 1 ) + ν(n M 2 − nM 2 ) (C5) It is integer in two limiting cases, confining holonomy ν = 1/3 at T < T c , and trivial holonomy, ν = 0 at large T . The unrestricted topological susceptibility is then