Confinement, Holonomy and Correlated Instanton-Dyon Ensemble I: SU(2) Yang-Mills Theory

The mechanism of confinement in Yang-Mills theories remains a challenge to our understanding of nonperturbative gauge dynamics. While it is widely perceived that confinement may arise from chromo-magnetically charged gauge configurations with nontrivial topology, it is not clear what types of configurations could do that and how, in pure Yang-Mills and QCD-like (non-supersymmetric) theories. Recently a promising approach has emerged, based on statistical ensembles of dyons/anti-dyons that are constituents of instanton/anti-instanton solutions with nontrivial holonomy where the holonomy plays a vital role as an effective"Higgsing"mechanism. We report a thorough numerical investigation of the confinement dynamics in SU(2) Yang-Mills theory by constructing such a statistical ensemble of correlated instanton-dyons.


I. INTRODUCTION
The Quantum Chromodynamics, or QCD, is established as the fundamental quantum field theory of strong nuclear force underlying all of nuclear physics. Despite its great success in describing an impressive variety of nuclear phenomena in nature, a key aspect of QCD remains mysterious and poses a great challenge to our understanding. While the theory has quarks and gluons as its basic degrees of freedom in the Lagrangian, the colored quarks and gluons are absent from the observed physical spectrum in which the various color-singlet hadronic states emerge instead. This phenomenon, often referred to with the broad term "confinement" (-see recent review in e.g. [1]), occurs also in a wide variety of QCDlike theories and notably in pure Yang-Mills theories. The latter fact makes it obvious that confinement arises from the nonperturbative gauge dynamics in the gluonic sector. It was suggested long ago [2][3][4], based on analogy with superconductivity, that the confinement may arise from chromo-magnetically charged and topologically nontrivial gauge configurations, with the vacuum being a "dual superconductor" of such magnetic objects (-see review in e.g. [5,6] ). This scenario is highly appealing and widely perceived to be a likely mechanism for confinement. The idea was extensively explored via lattice simulations as well as has been concretely shown to work in certain supersymmetric theories [7,8]. More recently, the idea of "dual superconductor" vacuum has been further advanced into a "magnetic scenario" for the hot plasma phase in the temperature regime above but near the transition temperature T c [9]: indeed, if there is a magnetic superconductor below T c , there should be a "precursor", i.e., the normal phase of thermal magnetic plasma just above T c . There are strong evidences from lattice simulations [10][11][12] for this scenario and most interestingly such a magnetic component is found to be crucial for explaining a number of key transport properties of the near-T c QCD plasma as measured from heavy ion collision experiments [13][14][15] (-see e.g. [16] for a recent review). Despite the various progress so far along this line of thought, it is nevertheless still unclear what types of configurations could drive confinement and how, in pure Yang-Mills and QCD-like (non-supersymmetric) theories. It shall be mentioned in passing that there exist a variety of interesting alternative ideas about possible mechanism and possible topological objects that may drive the confinement [17][18][19][20][21].
Let us elaborate a bit on the difficulty to identify the relevant topological configurations for confinement in these theories. Conventional instantons (as well as their finite-temperature counterpart, the calorons) [22,23] are well known topological objects and studied in detail (see reviews in e.g. [24,25]). However, these conventional instantons/calorons only have trivial holonomy, that is, trivial Polyakov loop at spatial infinity | x| → ∞ (-see Appendices A and B for a detailed discussion about holonomy), which is in sharp contrast to the confining vacuum where the holonomy is maximally nontrivial. Furthermore they are color neutral with no chromo-magnetic charge. Therefore, such conventional instantons/calorons can not be responsible for confinement. Another natural candidate would be 't Hoof-Polyakov type of magnetic monopole (-see e.g. review in [26]). However, a crucial difference of pure Yang-Mills or QCD (as compared with e.g. George-Glashow model or Seiberg-Witten therory) is that they do not have adjoint scalar fields which would provide the natural "Higgsing" on the spatial boundary, thus allowing the construction of magnetic monopole solutions. Because of this difficulty, one usually has to rely upon certain gaugefixing procedure to manifest the monopoles in the pure Yang-Mills or QCD cases.
As it turns out, both difficulties are resolved recently in a crucial new development: the construction of caloron solutions with nontrivial holonomy and nontrivial topology to the classical Yang-Mills equations, known as the KvBLL calorons [27][28][29] (-see Appendix C for a detailed discussion). First of all, such new calorons by construction acquire the necessary sensitivity to holonomy and thus are able to play a role in the confinement dynamics. Even more nontrivially, such a caloron is made of N c different constituent dyons (for SU (N c ) theory) which are chrome-magnetically charged and whose properties critically depend upon the holonomy. In these solutions, the nontrivial holonomy provides the needed nontrivial boundary constraints (which would have to come from the adjoint scalar fields in "Higgsed" theories). Owing to such important new features, the KvBLL calorons with their dyon constituents provide the unique topological configurations that could potentially account for the nonperturbative dynamics underlying confinement.
Based on the KvBLL solutions, a promising approach has emerged for understanding confinement in a statistical ensemble of dyons/anti-dyons arising from the constituents of the KvBLL calorons. Early works on an uncorrelated ensemble of these objects -to be called instanton-dyons (following recent literature) from now on in the present paper -already indicated that their contributions (alone) to the holonomy potential tends to push the system toward confining holonomy [30][31][32]. However it was later found that such contributions would not be enough to overcome the one-loop perturbative contributions to the holonomy potential (which favors the trivial holonomy -see detailed discussions in Appendix B). It was later found by Shuryak and collaborators [33][34][35][36][37] that an effective dyon-anti-dyon interaction with a short-range repulsive core (or in a broader term, a strong dyon-anti-dyon short-range correlation) appears necessary to enforce the confining holonomy at low temperature. Effective model for a Coulomb plasma of such dyons/anti-dyons was also constructed and shown to give a reasonable qualitative description of the low temperature T < T c properties of Yang-Mills theories [38,39].
In this paper, we report a thorough numerical investigation of the confinement dynamics in SU (2) Yang-Mills theory by constructing such a statistical ensemble of correlated instanton-dyons. We present high precision results for the temperature dependence of the holonomy potential, the order parameter for confinement transition, the dyon ensembe properties (densities and densitydensity correlations), as well as the temporal and spatial Wilson loops. In particular, we study the influence of the finite volume effect, the dyon-anti-dyon correlations as well as the screening mass on the confinement dynamics. Some of these results are considerably improved as compared with previous studies and many new results have not been previously reported.
The rest of the paper is organized as follows. In Section 2, we present the detailed construction of the correlated instanton-dyon ensemble as well as the numerical proce-dures. Our main results about the confinement dynamics in such an ensemble are reported in Section 3. The Section 4 then focuses on examining the consequences of the key parameters in the ensemble construction. Finally, we conclude the study in Section 5. In addition a few Appendices are included to explain some "background" information in detail and to make the paper more selfcontained for the convenience of readers.

A. The Partition Function
The construction of the partition function of the dyon ensemble begins with rewriting the one-loop quantum weight of a single KvBLL caloron in the limit of large dyon separation Eq. (D7) (-see Appendices C and D for more details), namely the contribution of a pair of L and M dyons, in the following way: Here, the fugacities per dyon species are introduced as f M = ΓS 2 e −νS ν 8ν 3 −1 and f L = ΓS 2 e −νSν 8ν 3 −1 . The P (ν) is the famous Gross-Pisarski-Yaffe result for perturbative contributions while the detĜ is the contribution of one-loop quantum fluctuation around one KvBLL caloron. Extending this result to arbitrary number of L and M dyons requires the inclusion of the appropriate metric of the moduli space which, as of today, its explicit form has not yet been found. Nevertheless, Diakonov and Petrov [30] proposed an approximate metric by merging that of a neutral cluster of dyons of different kind, namely, an L-M pair, with that of dyons of the same kind (originally proposed in [40]). Therefore, the full measure is approximated by the square of the determinant of a symmetric matrix G like det(g) ≈ det(G). Despite not being an exact solution, it possesses the interesting property that in the limit of K well separated L-M pairs, the measure factorizes into det(G) = det(Ĝ) K , i.e., as the product of K individual KvBLL caloron measures. It is thus straightforward to see that for a single L-M pair, the G matrix reduces toĜ (Eq. (D5)). In the SU (2) case, when the number of L and M dyons are N L and N M respectively, the dimension of this matrix G is (N L + N M ) × (N L + N M ) and its components are given by: where r mi is the position vector of the i th dyon of kind m (either L or M ). Furthermore, it should be clear that similar results can be obtained for anti-dyons.
As pointed out in [31], dyon-anti-dyon configurations are not saddle points of the Yang-Mills action. The inclusion of anti-self-dual fields in the ensemble is done by factorizing the measure into uncorrelated parts det(G D ) det(GD) (D for dyons andD for antidyons) times a correlated contribution e −V DD , where V DD is the action corresponding to dyon-anti-dyon interactions. Classical interactions between dyon-anti-dyon of the same kind was recently introduced in a gradient flow study in [35]. Using the parametrization found in [36], the potential takes the following form for ζ j > ζ c j and r jj = | r j − rj | . Below the limit ζ j < ζ c j , the interaction is repulsive and the proposed core potential for this region is given by where V c and ζ c j are the two key parameters that quantify the strength and range of the repulsive correlations between dyon-anti-dyon pairs.
Other interactions that have to be accounted for include the long-range forces due to the Abelian electric and magnetic charges and the nonlinear terms in the field strength tensor F µν , given by where e j and m j are the electric and magnetic charges (-see Table C.1) and h j = 1 for the M -type (anti)dyons while h j = −1 for the L-type ones. As expected, this gives exactly cancelled classical interaction between the L and M dyons (as well asL andM ) that together make a KvBLL caloron (owing to their BPS nature). On the other hand, there is repulsive interaction for the LM and ML pairs, while attractive interaction for the LL and MM pairs. In the construction of the ensemble, one has to sum over different number of (anti)dyons and also take into account the many-body screening effect which is introduced by means of a Debye screening mass M D . In doing so, all Coulomb terms appearing in the partition function (including those in the G matrices) are modified into r −1 → r −1 e −MDr . Combining Eqs. (3) to (5), the contribution from the inter-particle interactions to the action in the partition function is given by: With all the above elements, one finally writes down the following form for the full partition function of the dyonanti-dyon ensemble: where the factorial terms are needed to avoid duplicate counting of identical configurations with given numbers of dyons and anti-dyons. By requiring neutrality condition, i.e., equal number of dyons and anti-dyons of the same kind, the above expression can be further simplified into: where V is the 3D volume available and is obtained after performing integrals over all dyon positions. Finally, using Stirling's approximation log N ! ≈ N log N − N + log √ 2πN and defining the dimensionless dyon densities as n D = N D /V T 3 , we rewrite Z as a sum of weights running over number of dyons as In this framework, there are three key parameters as theoretical inputs: the screening mass M D , as well as the strength parameter V C and range parameter ζ c j for the dyon-anti-dyon interaction potential. In principle these parameters could be constrained by comparing relevant observables from the dyon ensemble with lattice simulations. Such quantitative comparison will be the goal of a forthcoming study, while the present paper focuses on qualitative question of demonstrating how the confinement is driven to occur in the correlated dyon ensemble.

B. The Monte-Carlo Simulations
Let us now discuss the details of the Monte-Carlo simulations to be used for evaluating the dyon ensemble partition function. Different from the implementation in [34,36,37], in our simulation we used a flat geometry, namely a box with periodic boundary conditions which shall be a more "realistic" approach and a more direct way to compare the results with e.g. lattice simulations.
From Eq. (10), it can be seen that all explicit dependence on the temperature T can be absorbed by rescaling rT → r, V T 3 → V , M D /T → M D and the free energy F/T ≡ − log Z → F . Since this simplifies the calculations, all the simulations are done using such scaled dimensionless variables. In doing so, the temperature T superficially disappears from the explicit simulations. However, the temperature dependence implicitly affects the system properties through the running coupling constant in the caloron action S, which at one loop level is given by (-see Appendix D) where Λ is the scale parameter in the regularization. Therefore, by varying S as a parameter in the simulation, one is essentially varying the system temperature.
It is straightforward to convert S into T /Λ. To further put temperature in e.g. MeV unit, one would then have to make a physical choice for the value of Λ. For example, one may choose Λ such that the critical temperature T c matches the lattice obtained value for SU (2) Yang-Mills theory. Once T c is fixed, one can then measure temperature T in terms of T /T c (-note this is equivalent to specifying the ratio T c /Λ).
The computation of all the observables are performed through Monte Carlo simulations using the well known Metropolis-Hastings algorithm. Each configuration is generated by first randomly varying the 3D positions of a single dyon or anti-dyon of each kind (and accounting for the periodic boundary conditions), then applying the acceptance algorithm, and moving to the next set of dyons/anti-dyons by repeating the same procedure. Once all positions have been swept, we then move to compute the observables with this new configuration and repeat all over again until the ensemble has been thermalized with enough statistics. It has been found that after about 2000 Monte Carlo configurations, the system is typically stabilized, after which the ensemble average would be calculated with the 10000 subsequent new configurations. The determined autocorrelation time was close to 5 configurations; therefore, the observables are averaged over 2000 configurations. On Fig. 3 we compare the free energy density calculated with a smaller number of Monte Carlo configurations for both confined and deconfined phases. The results are obviously consistent with each other and the small discrepancy between the two data sets is merely at a level of approximately 0.54% in the order parameter calculation L ∞ . This comparison clearly demonstrates that with 12000 configurations one obtains very reliable results with rather small statistical uncertainty.
One technical issue in the Monte Carlo sampling process is about the measure factors det(G). As it is an approximation to the actual moduli space metric, it may happen that some of the eigenvalues of the G matrices become negative thus violating the positive definiteness of the metric. This issue has been addressed in detail by [31,38,41]. To avoid configurations with negative eigenvalues in the simulations, the approach taken was to use the Metropolis-Hastings algorithm to reject such "wrong" configurations by assigning them a small statistical weight, i.e. if either G D or GD has at least one negative eigenvalue, then the weight exp(log det(G)) is substituted by e −100 , which was found to be enough to suppress these and to ensure an ensemble of sufficient configurations with all positive eigenvalues. It may be noted that this procedure effectively introduces a modification of the action, the impact of which is currently not well controlled and requires further investigation in the future.
One of the most important quantities to be calculated from the simulations is the holonomy potential or free energy density F/V at a given temperature. Due to the way it is defined, the calculation through Monte Carlo is not straightforward. However, there is a common method to evaluate it [35] which we will adopt here. Note that the only term that needs to be evaluated from the Monte Carlo configurations is e −V F since it is the only one depending upon spatial positions of the dyons/anti-dyons, while all other terms do not have such dependence. In the calculation, according to the standard thermodynamic integration, one introduces an auxiliary parameter λ as where and Dr is just the integration measure over all dyons' and anti-dyons' positions (for a total of 2(N L + N M ) of these particles in the simulation). It should be emphasized that for λ = 1, the above Eq. (12) is exactly equal to Eq. (9). The normalization factor V 2(NL+NM ) in the denominator above, is not introduced arbitrarily but rather comes directly from the construction of the partition function by correctly counting the "1/V" factors in the Eq. (9). This proper normalization factor also automatically gives F λ (0) = 0. Then, via standard Monte Carlo simulation procedure, one can compute the ensemble average of the following quantity: Lastly, by integrating out the λ dependence of the above, one arrives at the desired free energy: given that F λ (0) = 0 by definition. We emphasize again that for λ = 1, Eq. (12) reduces to Eq. (9), where the denominator V 2(NL+NM ) appears naturally from the construction of the partition function (Eqs. (7) to (9)), allowing to set F λ (0) = 0 unambiguously, regardless of the dyon numbers N L and N M . Finally, we discuss the choice of the parameters in this framework. For most of the results to be presented, we use a (dimensionless) spatial volume of the box to be V = 43.37. After several tests, it was determined that the optimal range of (anti)dyon density of each kind was n D ∈ [0, 0.5] (corresponding to N D ∈ [0, 22] number of (anti)dyons). Configurations with larger n D were found to have a rather small contribution to the partition function; therefore, discarded in the simulations (-see Section IV A). We choose the three key parameters as Debye mass M D = 2, the core potential strength V c = 20 and size ζ c j = 2; however, in Sections IV A to IV C, we will vary these quantities to explore the finite volume effects as well as the influence of the three key parameters.

A. The Holonomy Potential
It is known from lattice simulations that the SU (2) pure gauge theory has a certain critical temperature T c , with a confined phase at T < T c , a deconfined phase at T > T c , and a 2nd order phase transition connecting the two phases at T = T c . The relevant order parameter is the expectation value of the Polyakov loop at spatial infinity L ∞ (-see Appendix B) which is related to holonomy parameter ν by L ∞ = cos(πν). An expectation value of L ∞ = 0 or ν = 1/2 would correspond to the low temperature Z 2 center-symmetric, confined phase. A first important check is to examine whether such expected phase transition indeed occurs in the dyon ensemble. In order to see that, one needs to compute the holonomy potential, that is, the free energy density as a function of the holonomy F (ν) = −T log Z at varied temperature. Such holonomy potential determines the Polyakov loop dynamics and is a crucial input for a class of chiral models to incorporate confinement dynamics [42][43][44][45]. For any given temperature, the minimum of the holonomy potential determines the thermodynamically realized expectation value of the holonomy value which as order parameter thus tells us about the different phases of the theory. As mentioned earlier, the temperature dependence of all the observables in the ensemble comes from the instanton action S (Eq. (11)) which is an input parameter in the simulation. Fig. 1 shows the free energy density for S = 5, 6, . . . , 13. It is found that for S = 5 ∼ 7, the minimum of the free energy density lies at ν min = 0.5, namely maximal non-trivial holonomy corresponding to the confined phase. For S > 7, the shape of F/V becomes that of a symmetric double well potential with two minima located at ν min < 0.5 and ν min > 0.5 in a symmetric way. It shall be mentioned that as expected for an SU (2) pure gauge theory, F/V is symmetric under the interchange ν →ν = 1 − ν, and this feature has indeed been validated explicitly in the numerical calculations. So the results clearly reveal a confined phase at low temperature while a deconfined phase at high temperature.
To more accurately locate the critical action (or equivalently the critical temperature T c ), we further run the simulation for S = 7.25, 7.5 and 7.75. As shown on Fig. 2, for S ≥ 7.5 the Z 2 symmetry is clearly broken and the minimum of the free energy density is shifted away from the ν =ν = 0.5. For S = 7.25, more points were necessary to examine the minimum, and despite the potential on Fig. 2 exhibits a very flat dependence around ν = 0.5, the minimum was actually found around ν min ≈ 0.453. Thus, at the present numerical precision, we determine the critical temperature at S c = 7.22, which fixes our scale parameter from Eq. (11) at Λ = 0.373T c and allows us to express all temperature dependent quantities in terms of the ratio T /T c . We next come to the expectation value L ∞ = cos(πν), which can be determined from the position of the minimum of the holonomy potential. This is done by fitting the free energy density near the minima to a quadratic function with 9 to 15 points and then, through a derivative test on the fit, finding its minimum accurately.
As an important and insightful check of the role of L ∞ as an order parameter for the expected 2nd order phase transition, we quantitatively examine whether its dependence on temperature near the transition point follows the proper universality class. The well known Svetitsky-Yaffe conjecture [46], relates SU (2) pure gauge theory in (3 + 1) dimensions to the 3D Ising model of ferromagnetism by categorizing both in the same universality class, which has been proven several times in different numerical studies such as [47][48][49][50][51][52]. In this sense, L ∞ be- comes the analog of the magnetization, thus its critical behavior must follow the same universal power law with b and d the fitting parameters.
Using the well established values of the critical exponents of the 3D Ising model β ≈ 0.3265(3) and ∆ ≈ 0.530 (16) [53], on Fig. 4 we show the fitted curve obtained from the numerical results of the dyon ensemble in the near-T c region, namely 1 ≤ T /T c ≤ 1.274. The very low value of χ 2 = 1.44 × 10 −4 of the fit (which is partly due to the sizable error bar because of limited statistics) suggests an almost perfect agreement between the confinement/deconifnemnt phase transition behavior with the anticipated critical behavior of the 3D Ising model's 2nd order phase transition. It also demonstrates qualitative agreement with the lattice results from [54,55]. For completeness and comparison, we also show the fit using the mean-field critical exponent β mf = 1/2, which shows a qualitatively similar trend but a significantly larger value of χ 2 = 0.13. The comparison favors the former fitting and implies that the transition from the dyon ensembles captures the beyond-mean-field critical behavior of a 2nd order phase transition.
We now report the results for the expectation values of dyon densities, shown in Fig. 5. One can see that at T < T c , the L and M type densities are equal as expected. In the confined phase, the preferred holonomy corresponds to the maximally non-trivial one where both dyon types have the same core radius as well as equal action share and therefore equal weight in the partition  function. For T > T c , the prefered holonomy starts to shift away from the symmetric point towards the trivial holonomy (ν → 0 in this case) and the M dyons become larger and larger. Recalling from the KvBLL caloron solution (-see Appendix C), in the limit of trivial holonomy, the L dyon disappears and the field becomes that of the Harrington-Shepard caloron. A similar situation is observed in the ensemble as temperature is increased, with the L dyon density decreasing much faster than M type. The total density of all these magnetically charged objects demonstrates a strong temperature dependence with very rapid increase from high temperature toward near T c regime, in consistency with the magnetic scenario.

B. The Dyon and Anti-dyon Spatial Correlations
The interactions between dyons and anti-dyons are essential for the properties of the ensemble and in particular for driving the system toward confinement at high dyon density (i.e. low temperature) regime. The effect of such interactions can be illustrated by examining the spatial density density correlations among various pairs of dyons/anti-dyons, as defined in the following: which is normalized to that of an uncorrelated ideal gas and where the step function Θ δx (ξ) = 1 for 0 < ξ < δx and 0 otherwise. A value of unity for the G DD ′ would indicate a situation without correlations as is characteristic for a free gas ensemble. The numerical results are shown in Figs. 6 and 7 for all different (anti)dyon pairs combinations, each computed at one temperature value below T c and another one above T c . These results are obtained under equilibrium conditions for given temperature, i.e. with holonomy parameter ν being the one at minimal free energy and the number of dyons N D fixed at the ensemble averaged values. The presence of the repulsive core is clearly observed for all the dyon pairs, besides the LM for which there is none. At distances right above the core size ζ c j /2πν j the correlation functions seem to have a small bump that rapidly goes to unity at larger distance, indicating at a short-range correlation pattern arising from the repulsive core.

C. The Polyakov Loop Correlator
Besides the Polyakov loop expectation value itself, another important "indicator" of the confinement/deconfinement transition is the static (quark-antiquark) potential which essentially is evaluated from a temporal Wilson loop or equivalently the spatial correlator of the Polyakov loop. In particular the so-obtained static potential is expected, at large spatial separation, to exhibit a linearly rising behavior in the confined phase while to level off in the deconfined phase. It is important to evaluate this observable in the dyon ensemble.
The computation is however technically tricky in the present framework. In the large distance limit (| x| → ∞), the A 4 component of the dyon fields becomes Abelian (see Appendix C). However, the total A 4 of the ensemble (far away from their individual cores) cannot be given by a superposition of the individual fields of all dyons yet. Since the asymptotic condition A 4 | | x|→∞ = πντ 3 must be satisfied, one has to eliminate the holonomy parameter term in the gauge field associated with individual dyon by means of the time dependent gauge transformation U = exp(−iπνx 4 τ 3 ), after which one can then superimpose all dyonic fields and finally restore the asymptotic term with the inverse gauge transformation U † [35]. This procedure leads to where l( x) is the sum of all Coulomb terms of dyons and anti-dyons At finite temperature, the color averaged heavy quarkantiquark free energy F avg qq is defined through the expectation value of traced Polyakov loop correlators (-see Appendix B). For quarks in the fundamental representation, from Eqs. (18) and (19) and the definition of the Polyakov loop, it is straightforward to see that Thus, the color averaged static quark-anti-quark potential in the dyon ensemble is given by The above static potential, though, is different from a color-singlet static potential which is the one relevant for linear behavior at large separation. According to the color decomposition 2 ⊗2 = 1 ⊕ 3, an SU (2) quarkantiquark pair can interact through a singlet and a triplet channel [56], meaning that F avg qq is decomposed into where the singlet free energy is obtained from the following: and the triplet contribution follows trivially from Eq. (22).  Due to the periodic boundary conditions imposed in our geometry, the maximum allowed distance is | x − y| ≤ R/2, where R ≈ 3.51 is the size of the box of volume V = 43.37. To compute these observables, a total of 3000 Monte Carlo configurations are used for each combination of number of dyons (N L , N M = 0, ..., 22). To account for isotropy, for each interquark separation we averaged the contribution to the Polyakov loop correlator from 13 different orientations. At each temperature, the holonomy parameter ν is fixed to be the equilibrium value that the one which minimizes the ensemble free energy. In Figs. 8 and 9, we show the color averaged potential and its singlet and triplet contributions for the confined and deconfined phases at T /T c = 0.970 and T /T c = 1.674 respectively. In Fig. 10 we show the singlet channel free energy alone as a function of interquark separation | x− y| for several temperatures below and above T c . It may be noted that the color-averaged static potential above T c appears not fully saturated at large distance, due to two factors. The first is the finite volume effect (as will be discussed later in Section IV.A) which would limit the largest possible distance we could explore. The second is that in the high temperature deconfined phase, the perturbative thermal gluons (which are absent in the current framework) would contribute more and more importantly with increasing temperature to the screening of the static potentials.
An interesting comparison is with the static potential of quarks and anti-quarks in the adjoint representation, in which case no linear rising at large distance is expected as the gluons (themselves being adjoint) can screen out the potential. Following [44], one can obtain the adjoint static potential via the following relation with the fundamental one L f : where τ i are the Pauli matrices. Given that L ∈ SU (2),   the fundamental representation is generally defined as with a µ a µ = 1. Thus Eq. (24) can be rewritten as and it is easy to see that its trace is expressed in terms of the fundamental one as Therefore the adjoint static quark-antiquark free energy is then given by Notice we have included a normalization factor |Tr L a (0)| 2 in the correlator such that F a qq = 0 at | x − y| = 0. The resulting potentials are shown in Fig. 11 for different temperatures.
As pointed out already, at large separation in the confined phase, one expects the (fundamental representation) singlet static potential to have a linear rising behavior of the following form: where σ is the so called string tension. This is clearly observed in the fundamental representation (Fig. 10) at T /T c < 1. However, at temperatures above T c , the slope σ drops toward zero, as expected. In Fig. 12, we show the extracted string tensions for the singlet potential for several temperatures. In extracting the slope, we use linear fits for the large distance part but ignore the "curved" tails observed at the largest distances, which are most likely due to finite volume effects (-see Section IV A for an extended discussion). It may also be mentioned that our results show a relatively slow decrease of σ above T c , an effect which may also be due to finite volume issues. In contrast, the adjoint static potential does not show any linear rising at large separation. In short, our results from the dyon ensemble for the static quark-antiquark potentials in both representations are consistent with the expected behavior of an SU (2) pure gauge theory.

D. The Spatial Wilson Loop
Another interesting quantity to explore is the spatial Wilson loop. It is known that the spatial Wilson loop at finite temperature shows area law behavior with a finite spatial string tension σ s both below and above T c and thus by itself does not serve as an "indicator" of confinement transition [57][58][59][60]. Nonetheless, the restoration of Lorentz symmetry (Euclidean O(4)) at T → 0 suggests that in this limit, σ s should coincide with the string tension of the static potential Eq. (29) extracted from Polyakov loop correlators. The SU (2) traced spatial Wilson loop is defined as In the gauge where A 4 (x) is diagonal (-see Appendix C), the only non-vanishing spatial component of the dyon fields, in the asymptotic limit, is where m j = ±1 is the corresponding magnetic charge (Table C.1) and r = √ x i x i . The Dirac string singularity along the negative x 3 -axis, although a gauge artifact, might be an inconvenience for the numerical simulations. Therefore, for computing W C it is more suitable to use the corresponding magnetic field (instead of the gauge potential). For this, the Abelian Stokes theorem can be used to rewrite the spatial Wilson Loop in the so called "Abelian dominance" approximation [61] where B i ≡ 1 2 ε ijk F jk and A C is the area enclosed by a rectangular contour C. Therefore, the corresponding magnetic field to Eq. (31) is and B φ = B θ = 0. The total field strength from the whole ensemble will thus be It is interesting to examine whether the spatial Wilson loop computed from the dyon ensemble will follow the area law in both confining and deconfined phases, i.e.
In Fig. 13, the negative logarithm of W C in the fundamental representation is plotted as a function of the spatial loop area A C , which indeed demonstrates an almost linear rising behavior at large contour areas.
Recalling the units used in this work, the string tension obtained here is dimensionless after rescaling all quantities by temperature. To restore physical units, one makes the change σ s → σ s /T 2 . As has been established before [59], σ s increases with increasing T , however, σ s /T 2 should decrease with increasing temperature. Such a trend is consistent with our results from dyon ensemble. Finally, we've also examined the spatial Wilson loop for the adjoint representation, shown in Fig. 14. It is observed that the curve rises rapidly with loop area and reaches a plateau much faster than that of F a qq , again an indication of the screening effects for adjoint sources.

A. Finite volume effects
A rigorous study of all thermodynamic quantities in principle requires an infinite volume limit, which is obviously impossible for any realistic numerical simulations. In the case of the present study on the dyon ensemble, using a larger volume requires an increased number of dyons/anti-dyons in the simulations thus costing significantly more computing power. A practical approach would be to examine the finite size effect by perform tests with increasing volume of the box. In this Subsection we compare results obtained with two and three times the originally used volume, denoted as V 0 = 43.37.
One important feature to check is the (relative) contribution from various terms Z LM in the partition function Z Eq. (10). Note that for different volumes, each term with fixed number of dyons/anti-dyons N L,M would have different density. The better way to compare results computed with different volume would be to examine the contribution from given dyon/anti-dyon densities. In Figs. 15 and 16, we show how the contribution to partition function (from individual fixed-density terms in the ensemble sum) is changed with the increased volume. As expected, the maximum peak of the distribution becomes sharper around the most probable densities. Most importantly the location of the maximum does not change much with increased volume. Table I summarizes and compares numerical values of ensemble averages of the dyon densities at different temperatures as well as the free energy density for ν min obtained at the three different volumes. It can be seen that going from V 0 up to 3V 0 , there is a small shift (at few percent level) of the free energy density while small changes in the dyon densities. Such comparison clearly demonstrates that our thermodynamic results are quite stable with increasing system volume, which is an indication that our results shall be a very good approximation to the thermodynamic limit.
The finite volume also bears influence on the evalua-  (29) tion of spatial correlation observables, in particular the static quark-anti-quark potential. As mentioned in the previous section, it exhibits unnatural behavior near the largest distances that are allowed by the finite volume.
To test if this could indeed be a consequence of the finite volume, we've computed these observables with enlarged volume for comparison. In Fig. 17 we show the results of the singlet channel potential calculated in a box twice the volume of the original volume V 0 . For comparison, we also include the results from the original volume. One can see that indeed, the curved tails only appear at the edge of the box and at intermediate distances both potentials match substantially well. This comparison justifies our previous extraction of string tension via linear fit in intermediate distance regime and do suggest that for such spatial correlations, a significantly larger volume may be needed for their accurate evaluation.

B. The Influence of Dyon-anti-dyon Short-Range Correlation
A key ingredient in the confinement mechanism of dyon ensemble is the repulsive core potential V C jj . As defined in Eq. (4), there are two parameters which quantify such interaction: V c is the strength and ζ c j the size of the core. It is important to understand the influence of these parameters on the various observables. In Figs. 18 and 19 we show the free energy density as a function of ν for different values of V c at both low and high temperatures. A general observation is that a larger core strength V c always favors more the confining holonomy ν =ν = 1/2. A smaller V c , on the other hand, weakens the correlation and makes confinement harder to occur. Indeed for the V c = 10 case, even with the lowest temperature we explore, the system is still in the deconfined phase. These results also imply that the critical action S c needed for the confinement transition will shift toward larger values with increasing V c .
In a similar fashion, a change in the core size parameter ζ c j , will also result in considerable effect on the behavior As can be seen, when the core size is decreased, the free energy density's minimum shifts further and further away from the confining holonomy value ν =ν = 1/2. If the core size is too small then the system would be in deconfined phase even at the low temperature value computed here. With a large core size, the system could maintain a holonomy value near the confining one even at high temperature. The comparison clearly demonstrates the importance of the repulsive core. It is a strong repulsive core that drives the system toward favoring the confining holonomy at low temperature.

C. The Debye Screening Mass
Finally, we investigate another important parameter for the ensemble, namely, the Debye mass M D used to regularize the large distance behavior of the Coulomb terms and therefore to account for the screening effect. This parameter plays an important role in controlling the contributions to the free energy from the long range Coulomb interactions among the dyons/anti-dyons. To see its effect, we compare the free energy density versus holonomy from dyon ensembles with three different choices of the M D in Figs. 22 and 23 at both low and high temperatures. The results show that a smaller screening mass would disfavor the confining holonomy while a larger screening mass would help strengthen the confinement. This could be understood as follows: with a large screening mass the contribution to the free energy from many-body long-range Coulomb interactions get suppressed and thus the contribution from the short range correlations via the repulsive core, which essentially drives confinement, become relatively more important.

V. CONCLUSION
Confinement is a remarkable nonperturbative phenomenon in pure Yang-Mills and QCD-like theories. The mechanism of confinement remains a significant challenge to our understanding and is generally believed to be a consequence of certain nontrivial topological configurations of the gluonic sector. The recently found KvBLL caloron solutions with nontrivial holonomy, consisting of constituent dyons/anti-dyons, have provided a concrete and promising path of investigation. In this paper, we have constructed a statistical ensemble of such correlated instanton-dyons and performed a thorough nu- merical study of its various properties for the SU (2) Yang-Mills theory. Our main conclusion is that such an ensemble correctly produces the various essential features of the confinement dynamics from above to below the transition temperature. These features include the evolution of holonomy potential with temperature, a second order phase transition in terms of the order parameter (Polyakov loop expectation value), the linear static quark-anti-quark potential at large distance, etc. We have also found that the confinement dynamics is very sensitive to both the implemented short-range dyonanti-dyon correlations and the Debye screening effect in the many-body ensemble, by quantitatively investigating how the holonomy potential changes with these parameters. Given such success, it appears reasonable to believe that the ensemble of correlated instaton-dyons may indeed hold the key of confinement mechanism. The natural next steps of investigation would be the extension of the present framework toward the SU (3) case as well as toward the inclusion of dynamic fermions thus allowing the study of nontrivial interplay between the confinement transition and the spontaneous chiral symmetry breaking, which we shall report elsewhere in the future.
The classification of fiber bundles is an interesting topic in the geometry. Holonomy group, which describes the vector parallel transportation around closed loops, is one of tools to characterize the connection structure of a bundle. In gauge field theory the Wilson loop plays the same role as holonomy for gauge connections In the imaginary-time formalism of finite temperature field theory the temporal direction is compactified to a circle of radius (2πT ) −1 , where T is the temperature. Therefore, the holonomy could be defined around this loop as which is the so-called Polyakov loop. Here A 4 is an element of Lie algebra su(N ). And it could have different forms in different representations. Practically the Polyakov loop is very useful in the study of phase transition at finite temperature. The "condensate" of Polyakov loop serves as the order parameter of the confinementdeconfinement transition in the pure Yang-Mills theory.
To relate the fundamental quark confinement with the Polyakov loop it is intuitive to consider the free energy of a single static color charge, i.e., a quark with infinite mass [1] e −Fq/T = Tr q e −H/T /Tr e −H/T = n Ψ n (1 quark)|e −H/T |Ψ n (1 quark) n Ψ n (0 quark)|e −H/T |Ψ n (0 quark) .

(A3)
From the view point of path integral formalism the numerator is just an infinitely heavy quark propagating from ( x, 0) to ( x, 1/T ). Considering the kinetic part suppression due to the large mass the only contribution should come from the gauge field term which is equivalent to the Polyakov loop up to a constant e −Fq/T ∝ Tr L . (A4) Clearly, the quark confinement, which corresponds to F q = +∞, will induce Tr L = 0. Otherwise Tr L = 1 if the quark is totally free. In this sense the Polyakov loop could be treated as a order parameter for the confinement deconfinement transition in the pure Yang-Mills theory. And it is also apparent that the intension of confinement is much more than the condensate of Polyakov loop. More information could be revealed by studying the topological details of it.
In the fundamental representation at spatial infinity up to a global transformation [30,31] L = diag(e 2πiµ1 , e 2πiµ2 , ..., e 2πiµN ). (A5) Generators of SU (N ) group are all traceless, so eigenvalues should satisfy With a global transformation we could order them as And in this paper the term holonomy will be used especially to call the set {µ m |m = 1, 2, ..., N }. The holonomy is said to be trivial if L is in the center group Z N . Because Z N only has N one-dimensional complex irreducible representations, there are N choices for the trivial holonomy where k = 1, ..., N . Without triviality constraint there would be lots of choices for the holonomy. A typical nontrivial one is the so-called "maximally non-trivial" one Obviously, this is an equidistant one which yields to TrL = 0. Phase transitions often involve symmetry breaking or restoration. So it is for the confinement/deconfinement transition. With the definition of the Polyakov loop we can see that its condensate is the order parameter for the center symmetry. In finite temperature field theory the operation of the center symmetry is defined through the twisted gauge transformation, which satisfies the boundary condition along the imaginary-time dimension where z ∈ Z N which is in the center of gauge group SU (N ). It could be checked that the Yang-Mills lagrangian density is invariant under the center symmetry transformation. On the other hand, under such a gauge transformation U ( x, x 4 ), the Polyakov loop transforms as Using the boundary condition the trace of Polyakov loop should transform as In confined phase the Tr L = 0 means the center symmetry is preserved. While in deconfined phase it becomes nonzero which means the symmetry breaking.
Once the background configuration has been chosen we can always do the 1-loop perturbative calculation above this mean field. The essential part is to complete the integration and summation for dressed propagators of gauge fields [62,63].
It could be seen that the longitudinal part of δA i , which is an artifact of the gauge field, will cancel with the ghost part δA 4 [63]. Taking the group measure into account and performing the Matsubara summation, the holonomy dependent part of the perturbative potential energy is related to the integration and summation as where holonomy independent parts have been omitted at the second equation. And the Li 4 (z) and B 4 (z) are the polylogarithm function and Bernoulli polynomials of fourth order respectively. Here C represents different combinations µ j − µ k and β = T −1 . Gathering all of contributions from different combinations C = µ j − µ k , up to a holonomy independent constant the perturbative potential energy for the SU (N ) case is obtained as For the SU (2) case there is only one term ( The maximally non-trivial holonomy is obtained from the above equation, leading to where the Faulhaber's formula are used to complete the summation N i=1 i m . At this 1-loop level the perturbative potential energy has N minima corresponding to N elements of the center group. And the confining holonomy gives larger potential energy than the trivial ones. This means at 1-loop level trivial holonomies, which indicate the deconfinement, are favored at arbitrary temperatures. Hence, in order to achieve confinement at low temperature a more strict calculation is necessary with a topological non-trivial background configuration. The KvBLL Caloron is one of these choices.

Appendix C: The KvBLL Caloron Solution
The caloron field with non-trivial holonomy discovered by Kraan and van Baal [27,28] and independently by Lee and Lu [29] (therefore also known in the literature as the KvBLL caloron), is a classical solution to the SU (N ) Yang-Mills equations of motion in R 3 × S 1 . It is a selfdual field with unit topological charge and most importantly, the A 4 component can be gauged to be diagonal and constant at spatial infinity, which leads to a nontrivial Polyakov loop.
In the periodic gauge, the SU (2) KvBLL caloron field with period 1/T is given by Here,η a µν ≡ ε a muν − δ a µ δ ν4 + δ a ν δ µ4 are the so called 't Hooft symbols, T the temperature and τ a the Pauli matrices. The meaning of the s and r variables will be explained shortly. From this expression, it is not hard to see that at spatial infinity, the fourth component is indeed diagonal and constant A 4 | | x|→∞ = v τ 3 2 . This asymptotic value is parametrized as v ≡ 2πT ν, with ν ∈ [0, 1] and analogously,v = 2πTν withν = 1 − ν. Thus, the trace of the Polyakov loop at spatial infinity has the non-trivial form where ν = 1 2 corresponds to maximal non-trivial holonomy (L ∞ = 0) and ν = 0 trivial holonomy L ∞ = 1). Therefore, ν is naturally called the holonomy parameter.
The anti-self-dual caloron or anticaloronĀ µ with negative topological charge is easily obtained from Eq. (C1) byĀ As expected, the KvBLL reduces to the Harrington-Shepard caloron [23] in the limit of trivial holonomy (ν → 0 orν → 0). Furthermore, it becomes a standard BPST instanton [22] of size ρ in the zero temperature limit.
One of the most important properties of this solution becomes relevant when ρ ≫ 1/T . In this limit, the field is seen as composed of two constituent monopoles separated by a distance πρ 2 T . As ρ → ∞, the caloron becomes static and the monopoles are identified as the BPS type [64,65] with unit, but opposite, electric and magnetic charges, therefore named dyons or in this context intanton-dyons.
(Anti)dyons are commonly known as (anti)self-dual static solutions of the Yang-Mills equations of motion with an adjoint scalar (Higgs) field. However, one can construct dyonic solutions in pure Yang-Mills theory with the condition of non-trivial holonomy, namely A 4 | | x|→∞ = v. For SU (2) there are four kinds of dyon solutions which following the usual convention in the literature are labeled M and L for the self-dual fields andM andL for the antiself-dual ones, referred to as anti-dyons (-see Table C.1). In the hedgehog gauge, the M fields have the form of the common BPS monopole solution(for more details on the derivation refer to [31,66]) where n a = x a /| x| and the (lower)upper sign corresponds to the (anti)self-dual solution.
If a gauge configuration consists of more than two dyons, it is inconvenient to superimpose them in this gauge, since we are interested in configurations where all dyons have the same A 4 asymptotics at spatial infinity. This is achieved by using the matrices which satisfy the identity S ± (n a τ a )S † ± = ±τ 3 , and gauge-transform the dyon fields Eq. (C5) as In spherical coordinates, the dyon solutions in the new gauge take the form One should notice first that now the A 4 component is Abelian and equal for both M andM . Moreover, we have introduced a singularity along the negative x 3 -axis in A φ , a so called Dirac string which is merely a consequence of the gauge choice, hence the name stringy gauge.
The L andL solutions are obtained from Eq. (C8) by replacing v →v and apply two gauge transformations: first the time dependent U 1 = exp(−iπT x 4 τ 3 ) followed by a global rotation U 2 = exp(iπτ 2 /2) [29,67,68]. As required, these will leave the asymptotics of A 4 in the same form as for the M type solutions with the caveat that the spatial components are no longer static; however, in the large distance limit, neglecting exponentially small terms, the time dependent terms vanish and are no longer relevant in the scope of this article. The L type dyon fields in the stringy gauge thus are Going back to the KvBLL field, the emergence of such configurations suggests to express the caloron in terms of the "constituent" dyon's positions. The coordinates used to write the caloron in Eq. (C1) are then the positions of the dyon's center of mass denoted by r L and r M , the dyon separation r LM ≡ | r L − r M | = πρ 2 T , which for convenience is chosen to be along the x 3 -axis (-see This monopole picture is more evident when looking at the caloron in the vicinity of one of its constituent dyons and far away from the other, namely at large separations. For instance, near the L dyon center and far away from the M dyon (s ≫ 1/v), the caloron field reduces to that of the L dyon, whose asymptotic behavior is given by (see Eq. (C9)) In broad terms, this semiclassical procedure consists in taking the classical solution as a background field such that the gauge fields in the functional integral are where a(x) is a small quantum fluctuation of the classical solution (the KvBLL field). Then expand the action around the saddle point up to the desired order in a µ and compute the functional integral.
In [70] Diakonov et al. obtained an analytic expression for the quantum weight of the SU (2) KvBLL caloron in the one-loop approximation. They showed that in the limit of large separation between the constituent dyons (in the temperature scale) r LM ≫ 1/T , it can be written as Λ PV e γE 4πT 22 3 1 T r LM 5 3 (1 + 2πT ννr LM ) where P (ν) = (4π 2 /3)T 4 ν 2ν2 is the one-loop perturbative potential [62,63] (-see Appendix B), C ≈ 1.03142 is a combination of universal constants and the linear term in r LM proportional to P ′′ (ν) from the exponential factor has been ignored in this work. This expression can be further simplified in the approximation where the separation between dyons is much larger than their core sizes r LM ≫ 1 2πT ν , 1 2πTν ; taking the form Z KvBLL = e −V P (ν)/T d 3 r L d 3 r M T 6 (2π) To obtain Eq. (D3), one has to calculate the invariant measure of the moduli space metric of the caloron field denoted as det(g). In the general case of SU (N ), this is shown to be exactly equal to the determinant of a N × N matrixĜ [30,71], which for SU (2) is given by det(g) = det(Ĝ), which in the limit of large dyon separation reduces to det(Ĝ) ≈ 16π 2 νν, and thus the partition function Eq. (D4) is rewritten as Z KvBLL = e −V P (ν)/T d 3 r L d 3 r M T 6 (2π) where we have absorbed all constants into Λ. At the one loop calculation, the g −8 coupling in Eq. (D7) is not renormalized; however, a two loop improvement (ignoring the effects on P (ν)) will give