Chiral Symmetry Breaking and Confinement from an Interacting Ensemble of Instanton-dyons in Two-flavor Massless QCD

In this work we present the results from numerical simulations of an interacting ensemble of instanton-dyons in the $SU(3)$ gauge group with $N_f=2$ flavors of massless quarks. Dynamical quarks are included via the effective interactions induced by the fermionic determinant evaluated in the subspace of topological zero modes. The eigenvalue spectrum of the Dirac operator is studied at different volumes to extract the chiral condensate and eigenvalue gap, with both observables providing consistent values of the chiral transition temperature $T_c$. We find that a sufficient density of dyons is responsible for generating the confining potential and breaking the chiral symmetry, both of which are compatible with second-order transitions.


I. INTRODUCTION
A. Instantons-dyons at finite temperature Quantum Chromodynamics (QCD) possesses approximate symmetries whose breaking/restoration correspond to certain phase transitions. At finite temperature and zero chemical potentials, QCD has two crossovers: deconfinement and chiral symmetry restoration. In the infinite quark mass (pure gauge) limit, the theory has an exact Z 3 symmetry, for which the average Polyakov loop P is an order parameter. (The Polyakov loopP is a matrix and we define the average Polyakov loop as the scalar quantity P = 1 3 T r[P ( x)] .) At finite quark masses, it is no longer a strict order parameter, but still has a rapid change and a psuedocritical temperature that can be identified from its inflection point.
The average Polyakov loop is associated with the deconfinement of quark degrees of freedom by its relation to the heavy quark free energy [1] P = exp(−F Q /T ). (1) Of course, at large T the quarks are asymptotically free with trivial Polyakov loop P → 1. The question is then: which nonperturbative interactions drive the theory to the confining Polyakov loop P = 0 at low temperatures?
The instanton in SU (N ) gauge theories is the minimum of the action with vanishing fields A a µ at space-time infinity. At finite temperature, with a nonzero VEV of the Polyakov loop, one component along the Euclidean time direction is nonzero A 4 = 0. Looking for solutions of Yang-Mills equations with such modified conditions at infinity, it was found that the instanton dissolves into N c constituent dyons (also known as instanton-monopoles) connected by Dirac strings [2,3]. Unlike the instantons, the dyons interact directly with the holonomy. It was suggested then that the dyons can generate a confining potential which overcomes the perturbative interactions of thermal gluons. For a review, see e.g. Ref. [4].
Analytic descriptions of how the dyons generate confine are possible in particular supersymmetric models [5,6], which can be achieved with a dilute gas of dyons due to a cancellation of the deconfining potential. In the standard Yang-Mills theories, the deconfining Gross-Pisarski-Yaffe potential [7] means that a dense, stronglycoupled ensemble is required to confine the theory. In these cases, the dyons have been shown to generate confinement via numerical simulations in both the pure SU (2) [8,9] and, more recently, SU (3) [10] cases.
One can also take the 'inverse' approach, identifying dyons in lattice configurations using the Dirac operator and gradient flow methods to reveal the dyons [11][12][13] from underneath the quantum fluctuations of the gluon field.

B. Instanton-dyons and fermions
If two light quark flavors are massless, QCD has an exact SU (2) L × SU (2) R chiral flavor symmetry. Below some T c the axial part of this symmetry is spontaneously broken. The chiral quark condensate qq serves as the order parameter of the transition.
The interactions between instantons and quarks have been studied since the 1980's, starting with the instanton-induced 't Hooft Lagrangian. This interaction explicitly breaks the U A (1) symmetry via the topological quark zero modes. Later, the Instanton Liquid Model (ILM) [14] showed that the breaking of chiral symmetry is related to the collectivization of said zero modes into the so-called zero-mode zone (ZMZ). For a review, see e.g. Ref. [15].
Following the discovery of the instanton dyons, it was shown that the quark zero mode of the instanton localizes to a single constituent dyon. Which one depends on the quark periodicity condition [16]. These zero modes, like the dyons themselves, have explicit dependence on the holonomy. In the case of (physical) antiperiodic quarks, all zero modes are localized to the L dyons as will be discussed in the next sections.
Previous studies of ensembles of dyons have analyzed the Dirac eigenvalue spectrum in this zero-mode zone and confirmed their role in chiral symmetry breaking in the arXiv:2108.06353v1 [hep-ph] 13 Aug 2021 SU (2) gauge group via mean-field methods [17] and numerical studies [18,19]. Both techniques have also been employed to study the phase transitions in theories with modified quark periodicities [20,21]. The goal of this work is then to extend such numerical studies, in an effort to approach physical QCD, to the SU (3) gauge group with N f = 2 flavors of massless, dynamical quarks, looking at both the confinement and chiral symmetry transitions.
More recently dyons have been identified on the lattice via their quark zero modes [13]. It has been shown that the form of the quark zero modes are remarkably insensitive to the many perturbative gluons in which they are submerged. From studies such as this, dyon densities and correlation functions can be calculated, serving as useful constraints on models of the dyon interactions such as ours.
The structure of this paper is as follows: In Section II the physics of the dyon ensemble, their interactions, and the holonomy is described. Section III focuses on the quark zero modes and the interactions induced by them. Technical details of the simulations are discussed in Section IV. Finally, Sections V and VI lay out the results relevant to the deconfinement and chiral phase transitions, respectively.

A. Holonomy and the dyon partition function
The instanton-dyons are obtained by generalizing the instanton solution to nonzero holonomy (nontrivial Polyakov loop) [2,3]. In SU (3), the instanton is decomposed into three dyon species, the M 1 and M 2 dyons corresponding to the maximally diagonal subgroup and the 4-time dependent L dyon, as well as corresponding antidyons. The holonomies are the differences in the phases µ 1 , µ 2 , µ 3 of the eigenvalues of the gauge field component A 4 at spatial infinity. The holonomies are defined as ν i = µ i+1 − µ i with ν i = 1. Demanding that P be real reduces the individual dyon holonomies to depend on a single parameter ν.
The dyon actions and core sizes are directly related to their individual holonomies. The actions of each dyon, in terms of the instanton action S 0 , are The sizes of the dyon cores scale as 1/ν i . We work in dimensionless units with all lengths in units of 1/T , and define the following dimensionless quantities: the volumeṼ 3 = (LT ) 3 , free energyF = F/T , free energy density f =F /Ṽ 3 , and dyon densities n i = N i /Ṽ 3 . The single holonomy parameter is defined on the interval ν ∈ [0, 1/2] and is related to the Polyakov loop via The main scale of interest, the instanton action S 0 , is related to the temperature by where N c = 3, N f = 2, and Λ = 4.8. The value of Λ is chosen to set the critical temperature to be around S 0 ∼ 12. The dyon partition function and interactions are identical to that of the ensemble in the pure SU (3) theory [10] with two additional terms: a potential arising from the perturbative interactions of the quarks with the holonomy, and quark-induced interactions between L andL dyons. Here we provide a brief description of the partition function, focusing on the new terms.
The partition function of the dyon ensemble is separated into two parts as Z = Z 0 Z int , where Z 0 is the non-interacting terms and Z int contains all of the dyon interactions. The partition function is where d ν is the weight of an individual dyon with holonomy ν [22] The potentials in the partition function are the perturbative potentials of the quarks and gluons with the holonomy. The gluons experience the well-known Gross-Pisarski-Yaffe potential V GP Y [7], which in SU (3) is The first term new to this work is the second potential V quark , which is the analogous interaction of N f flavors of massless quarks with the holonomy [23,24]. For SU (3), it has the form Both of these potentials favor the deconfining holonomy ν = 0. Fig. 1 shows the contributions of both terms to the free energy density.

B. Dyon interactions
Within the temperature range of interest near the critical temperature T c , the dyons are strongly interacting. These interactions contribute significantly to the partition function and their computation is nontrivial, requiring Monte-Carlo integration over the 3N D -dimensional space of the dyons' collective coordinates. The interaction terms of the partition function can be written in the form Here the interactions are separated into three parts: the classical binary interactions of the dyons ∆S cl , the oneloop fluctuation determinants det (G) and det (Ḡ), and the eigenvalues λ i of the Dirac operator due to the inclusion of N f flavors of dynamical massless quarks. The classical and one-loop interactions included are the same as those in the pure SU (3) theory. At large distances, the dyon interactions are Coulomb-like. We use the parameterized form between dyons i and j. The coefficient C ij is −2 for dyonantidyon pairs of the same type, 1 for dyon-antidyon pairs of different types, and 0 for dyon-dyon or antidyonantidyon pairs. At distances shorter than the core size x 0 = 2πν i r 0 T , dyons of the same type, regardless of duality, experience a repulsive core of the form ∆S core cl = νV 0 1 + e 2πνT (r−r0) .
The volume metrics (G for the dyons andḠ for the antidyons) are the Diakonov determinants [4], each with the same form for the elements between the i-th dyon of type m and the j-th dyon of type n The other term in the interactions is fermionic determinant, not included in our previous work. It can be written as the product of the eigenvalues of the Dirac operator det( / D) = λ i . Its form is discussed in detail in the next section.
Additionally all Coulomb-like terms in the classical and one-loop interactions are regulated by an electric Debye screening term e M D rT . There are three phenomenological parameters of the model which at present are not known from first principles, V 0 , x 0 , and M D . We use the values V 0 = 40, x 0 = 1.8, and M D = 1.5 in this work.
The contribution of the interaction potential to the free energy density is computed via the standard integration over a dummy parameter λ, In total, the free energy density for an ensemble with specified input parameters is where Z 0 has been expanded with Stirling's approximation to three terms. The central goal of this work is then to compute f for a range of input parameters and determine the location of the minimum for each value of T , thereby determining the physical properties of the ensemble as functions solely of the temperature.

A. Fermionic determinant
The main result of incorporating dynamical quarks into the gauge theory is the inclusion of the fermionic determinant in the partition function. In the context of the instanton-dyon ensemble, the fermionic determinant is approximated by considering only the subspace of quark states spanned by the zero modes. This approximation results in the 'hopping matrix' assuming equal masses for all quark flavors a. Here we have removed a factor of i from the l.h.s. of the equation in order to writeT as a purely real matrix. For physical (antiperiodic) quarks, the right-handed zero modes are localized on the L dyons and the lefthanded zero modes are localized on theL dyons. Individual elements of the hopping matrix can be interpreted as the 'hopping amplitude' for a quark going between L dyon i andL dyon j. These elements are then given by (16) If one approximates the total gauge field as the sum of the fields of two dyons, then the covariant derivative can be reduced to an ordinary derivative by the zero mode equations of motion. This hopping amplitude only has a nonzero contribution from / D when the zero modes have opposite chirality (i.e. hopping within LL pairs, but not LL orLL) and because the zero modes are orthogonal, the mass term only contributes along the diagonal. Thus, the hopping matrix takes on the simple form The hopping matrixT is a 2N L × 2N L antisymmetric matrix in the case of massless quarks which we consider in most of this work. The fermionic determinant may be interpreted as a sum of all closed loops of the quarks hopping between dyon-antidyon pairs. As an explicit example, let us consider the case of two L dyons 1 and 2 and twoL dyons1 and2. In this case the determinant of the hopping matrix is The first two terms correspond to the two possible twoloop diagrams in which each L dyon forms a closed loop with a singleL dyon each (the upper diagrams in Fig. 2). The last term represents the one-loop diagram in which the quark hops between all four dyons (the bottom diagram in Fig. 2). For arbitrary number of L-dyons N L , there are N L ! terms in the determinant of the hopping matrix. (Of course, numerically we evaluate the determinant by a more efficient algorithm.) Whether the determinant is dominated by paths involving single pairs ('instanton molecules') or paths of many dyons ('instanton polymers') determines the state of the chiral symmetry. Consider the case in which the dyons are arranged into well-separated LL pairs. In this case, the determinant is dominated by T ii and T ij → 0 for i = j. In this case the eigenvalues of the matrix are λ ≈ ±T ii . In this configuration, the eigenvalues are then all large for the nearby dyon-antidyon pairs, suppressing the density of near-zero eigenvalues, leading to the disappearance of the quark condensate. It is in the collectivized 'polymer' regime, where there are many T ij 's of comparable magnitude which nearly cancel each other in the eigenvalues, which possesses a non-vanishing density of eigenvalues near zero.

B. Parameterization of the hopping matrix elements
The general form of the quark zero modes on the SU (N ) dyon gauge fields was first given in Ref. [16]. An explicit form for the zero mode density, in terms of all dyon coordinates is far too complicated to be of a practical use. (See Appendix A for a discussion of the general solution.) Instead, we start here with the form of the zero mode for a lone L dyon 1 : where the symbol contains the color and spin structure. The left-handed zero mode on theL dyon is found by changing the spin of the quark in the symbol. The spatial structure of the wavefunction is the same for all SU (N ). In this limit, the zero mode density has no Euclidean time dependence. The effect of nearby M i dyons is to destructively interfere with the gauge field of the L dyon localizing the zero mode in both space and time.
For some observables, such as hadronic correlation functions [25], this inon-binary forces seem to be necessary to achieve reasonable results. In some preliminary computations for this work, an adhoc approach to including the effects of interference was considered by appending a term estimating the localization effects to the wavefunction in Eq. (19). It was found that these terms led to only a modest modification of the quark-induced interactions. Thus we do not include such effects in out parameterization of the hopping matrix.
Computing T ij requires numerical integration, and instead a parameterization of it must be used. For a detailed derivation of T ij see Ref. [19]. We use the same parameterization as in the SU (2) work [18], known there as 'Parametrization A,' The magnetically-charged dyons have Dirac strings which are, in principle, pure gauge artifacts. However, the use of the sum ansatz in combining the two gauge fields introduces gauge-dependent factors in the zero modes. This leaves ambiguity in the overall normalization, handled here by c which we treat as a tunable parameter of the model and have chosen to use ln(c ) = 4.45.
The fermionic determinant adds an effective potential to the ensemble of the form This induces complicated many-body interactions between all L andL dyons, but at the simplest level is an attractive force within LL pairs.

IV. SIMULATIONS AND ANALYSIS
The simulation setup and analysis follows much of what was done in the previous work (see Section III of Ref. [10]). The free energy density of the dyon ensembles are computed via Monte-Carlo integration using the standard Metropolis algorithm. For each set of input parameters the simulation is run with 10 values of the dummy parameter λ = 0.1, ..., 1. Every dyon position is updated five times between samplings and 2000 configurations are sampled for each value of λ.
Each simulation is run at fixed dyon number N D = 120 (60 dyons and 60 antidyons). The densities of the dyons are controlled by varying the length L of the sides of the simulation box and the relative number of each type of dyon. Periodic (spatial) boundary conditions are imposed by a set of 26 image boxes placed around the main simulation box. Because of the large cost of computing the determinants of the matrices, only dyon interactions within a local box of length L centered on the dyon whose position is being updated are computed at each update step. Compared to the pure SU (3) case, the only addition to the Metropolis update step is computing the fermionic determinant. Because it only involves the positions of the L andL dyons, it is much smaller than the Diakonov determinant and does not increase the computational cost significantly. This is in contrast to lattice simulations, where the fermionic determinant is typically the most computationally-expensive task.
For any finite number of images there are dyons near the faces of the total simulation box which feel unphysical effects of the boundaries. In the pure SU (3) case, where all interactions were exponentially suppressed by the Debye mass, these finite-volume effects were only a percent-level correction to the free energy. The quarkinduced interactions are linear at long distances and must be treated more carefully. When computing the potential from these interactions, for each pair of dyons i and j one should consider only the image of j which minimizes the distance between the pair. This ensures that the short-distance (large-eigenvalue) interactions of dyons near the faces of the boxes aren't excluded while long-distance (small-eigenvalue) of dyons on the opposite sides of the boxes aren't over-counted. This introduces an effective cutoff distance L/2 considered in the hoping matrix which is remedied when taking the infinite-volume limit. Without using this technique, the eigenvalue distribution has a large number of unphysical, arbitrarilysmall eigenvalues. The physical parameters of the dyon ensemble are determined by a fit near the free energy minimum in the space of the input parameters for each value of S 0 . Once the minima are found, additional simulations are run with the physical parameters generating 60,000 configurations to compute the eigenvalue distributions with better statistics. In order to more accurately represent both fitted densities simultaneously, the dyon number of these runs is allowed to vary slightly with 120 ≤ N D ≤ 128. Additionally, in order to study finite-volume effects, 30,000 configurations are generated with 2N D dyons and 20,000 configurations are generated with 3N D dyons at the physical parameters for each temperature (thus the number of eigenvalues computed is the same for all three ensemble sizes).

V. THE POLYAKOV LOOP AND DECONFINEMENT
The physical properties of the dyon ensemble are determined by the location of the free energy minimum in the space of input parameters. The holonomy-and densitydependence of the potential determines the Polyakov loop and the nature of the deconfinement phase transition. Fig. 3 gives an example of the holonomy potential. The same general features are seen here as in the pure SU (3) case, namely that at high densities the minimum is located in the confining phase (ν = 1/3) as the dyon interactions dominate, while as the densities are reduced, the minima are pushed to lower values of ν, driven by the perturbative potentials. It is the location of the global minimum, at some intermediate densities, that represents the physical value of the holonomy in the thermodynamic limit. Crucially, however is the fact that the holonomy now varies smoothly as a function of temperature. Unlike the pure SU (3) case, there are not two nearly-degenerate minima near T c (see Fig. 4 of Ref. [10]). In QCD the deconfinement transition is a smooth crossover occurring almost simultaneously with the chiral transition, T deconf ∼ T c [26,27]. In the N f = 2 massless case we consider, evidence is less conclusive. Recently much progress has been made, in particular by the Bielefeld group [28][29][30][31][32][33], on studying the phase transition in (2 + 1)−QCD on the lattice in the chiral limit -QCD with massless up and down quarks and a physical strange quark m s ≈ 95 MeV. Of note is the fact that the location and form of the deconfinement transition is very sensitive to the light quark masses.
The order of the phase transition in the chiral limit is dependent on the state of the U (1) A symmetry breaking near the chiral restoration temperature. If the U (1) A breaking remains significant at these temperatures, as is expected, the phase transitions are expected to be second order belonging to the O(4) universality class, for large enough values of the strange quark mass m s [34,35]. If the U (1) A breaking is small near T c , the transition may be first order, although the recent lattice results disfavor that scenario. The phase diagram of the O(4) class is characterized by the temperature T and external field H. In the case of the gauge theory, the role of the external field is played by the light quark masses, meaning we look at the H = 0 line with the massless quarks. We then fit the data for P (T ) a the form inspired by Ref. [33].
where t = (T −T deconf )/T deconf and α = 2−β(1+δ) is the hyperscaling variable with the values of the O(4) class β = 0.380, δ = 4.824, and α = −0.2131 [36,37]. From this fit (χ 2 = 0.0257) we find the critical temperature of the deconfinement transition S 0 (T deconf ) = 10.44 ± 0.29. We also constrain the signs of the fit parameters to match those determined from the lattice data a 0 , A > 0, a 1 < 0. Removing these constraints or treating α as a fit parameter results in comparably-good fits with very different T deconf values. For example, the O(2) universality class is qualitatively similar but α = −0.0172 [38,39] is an order of magnitude different from the O(4) value. The dyon data shows slightly better agreement to that fit (χ 2 = 0.0178), so we do not claim that our data supports O(4) over other O(N ) universality classes, but merely that it is compatible with the expected O(4) behavior. Because of the multiple potential fits, one should consider the determination of T deconf to have larger uncertainties than those given by any individual fit.
In QCD, with no exact chiral or Z 3 symmetries, the transitions are an analytic crossover. The psuedocritical temperature is defined by the inflection point of the curve, where the derivative with respect to the temperature has a maximum, since there is no universal scaling expected near the transition from which a proper critical temperature could be determined via a fit. Lattice QCD studies find that, for the Polyakov loop, the peak is broad and hard to accurately define. We determine the derivative of P (T ) by simply computing the slopes between points. We see in Fig. 4 that no distinct peak can be effectively determined due to the large uncertainties.
The dyon densities, both shown in Fig. 5, decrease as the temperature rises. Unlike the pure SU (3) case, the L-dyon density remains smooth around T c . Even in the confined phase, when the dyon actions and sizes are equal, the densities are not due to the quark-induced interactions breaking the symmetry between the dyon types (or more generally, breaking Z 3 symmetry). The ratio of the densities at low T is directly sensitive to our choice of c .

A. Eigenvalue distribution of the Dirac operator
The breaking of chiral symmetry is associated with the existence of a nonzero quark condensate qq generated by nonperturbative effects at low temperatures. The quark condensate is related to the zero eigenvalues of the Dirac operator by the Banks-Casher relation [40] where ρ(λ) is the spectral density of the Dirac operator.
With the hopping matrix being antisymmetric, its spectrum is symmetric in the sign of the eigenvalue, ρ(λ) = ρ(−λ). For simplicity, we only show the positive spectrum in the plots. Fig. 6 shows examples of the spectrum in both phases.
The smallest eigenvalues correspond to collectivized modes with the zero modes of many dyons overlapping. For a system of finite size, the eigenvalues below some value λ min are suppressed with λ min ∝ 1/V , meaning that regardless of which phase the system is in, a finite system will always have zero eigenvalue density near zero.
In Fig. 6 (left), the steep decrease in eigenvalues below λ ∼ 0.06 is a result of the finite size of the system. The dashed line extrapolating to λ = 0 gives an estimate of the spectral density in the absence of such effects. In the restored phase (Fig. 6 (right)), the there is a larger range of decreasing eigenvalues. The decrease in density near zero reflects a real absence of collectivization rather than finite volume effects. This can't be seen from a single eigenvalue distribution alone, and distinguishing the two phases requires analyzing the distributions for multiple volumes.

B. Infinite-volume extrapolation and results
We locate the chiral phase transition by two different methods. The first is by extracting the chiral condensate Σ(T ) by extrapolating the small eigenvalue distribution to infinite volume using results from random-matrix theory. The vanishing of the condensate is related to the psuedocritical temperature T c . The other way is to use the eigenvalue gap ∆(T ) by fitting the smallest eigenvalues to a linear function (similar to Fig. 6 (right)). Above a temperature T gap , the restoration of chiral symmetry leads to a finite size of the smallest eigenvalues, meaning the lowest excitations have finite mass.
The mesoscopic volume scaling of the near-zero eigenvalues can be understood from random-matrix theory, which for N c = 3, N f = 2, gives a Dirac eigenvalue spectrum of the form [41]  where z = λV Σ 1 and J n are the Bessel functions. Here Σ 1 is the scaling factor and Σ 2 is the overall normalization factor. In the infinite volume limit ρ(0) → V Σ 2 .
Both of the factors Σ 1 and Σ 2 are related to different physics with different volume-dependence. In the case of the dense, low-temperate ensemble the eigenvalue distribution should be dominated by the collectivized modes. In this case, increasing the volume should reduce the region of suppressed eigenvalues by the same factor and Σ 1 ∝ V ∝ N D . On the opposite end, when the ensemble is dilute and comprised of dyon-antidyon pairs, Σ 1 is independent of the system size. Of course, in the region near T c , the ensemble is a mixture of both components and we must interpolate between the two.
A fit to the distribution gives two parameters per ensemble size. To extract the infinite-volume value of the condensate we use an interpolating function to determine how much of the region of smallest eigenvalue is real or a finite volume effect. We use the function where Σ 2 = (Σ V 2 + Σ 2V 2 + Σ 3V 2 )/3 is the overall scale and each term in parenthesis we call a scaling factor.
The philosophy of these scaling functions is as follows: each scaling factor is linear in the ratios of Σ i 1 and chosen such that it gives 0 or 1 in the opposite cases described above. In the case where the lowest-eigenvalue portion of the spectrum is real and doesn't change with the volume, the scaling factors are 0, and thus the condensate is 0 as there is a finite eigenvalue gap. In the case where the suppressed region is entirely due to finite volume effects, the factors should scale with V and Σ 3V 1 = 3/2Σ 2V 2 , Σ 2V 1 = 2Σ V 2 the scale factors are 1 meaning the overall value of the condensate remains in the infinite-volume limit Σ = Σ 2 . Additionally we enforce an upper bound of 1 on each scaling factor in cases where the values scale faster than V . Two scale factors are used to scale between the three volumes. Alternative functions using just two of the volumes are discussed in Appendix B.
Along the H = 0 line, the chiral condensate (analogous to the magnetization M in the O(4) spin model) takes the form where C is some (non-universal) constant. We qualitatively compare this form to our results in Fig. 8. As expected, we see a very rapid drop in the condensate to zero just below T c . More data points just above and below T c would be needed to get a more accurate fit to the data. Going to the lowest temperatures we do not see the continued increase in Σ expected by Eq. (26), however this universal scaling behavior is only applicable in the region near the critical temperature. As with the deconfinement transition, we claim our results are compatible with the second-order transition of the O(4) class. Certainly, like the lattice data, our results disfavor the potential first-order transition.
To determine the eigenvalue gap, a linear fit is performed on the smallest eigenvalues for each of the ensemble sizes. The x-intercept of each fit gives the eigenvalue gap. In order to extract the true value of the gap and remove finite-volume effects, a linear function in 1/V is fit to the three gaps for different volumes. The value of the gap at 1/V = 0 gives the true value of the gap ∆(T ) in the absence of the finite-volume suppression of small eigenvalues. The results of this analysis are also seen in Fig. 8.
Below a certain temperature, the gaps are all compatible with zero, indicating a finite density of zero eigenvalues. At S 0 ∼ 13.0 − 13.5, a finite gap forms and quickly rises, nearly simultaneously with the steep drop in the value of the condensate. Both of the fitting methods for the condensate Σ(T ) and gap ∆(T ) give consistent temperatures for the chiral symmetry restoration. From a fit of the condensate to the O(4) form, we get S 0 (T c ) = 13.14 ± 0.14, the central value of which sets the relative temperature scale in our plots. With this result, there is still one nonzero data point above T c . The fit shown in the plot provides better agreement with the data points just below T c and is also more in line with the appearance of a nonzero gap. Like deconfinement, in massive QCD the chiral phase transition is an analytic crossover [42], but is second order in the massless case. One of the most interesting results from the (2 + 1)−QCD lattice works [28][29][30][31][32][33] is that taking the chiral limit and reducing the light quark masses by just a few MeV results in a significantly-reduced tran-sition temperature T c [32,42]. Critical temperatures for the physical, massive case and the massless limit are T QCD c = 156.5 ± 1.5 MeV, T 2+1 c = 132 +3 −6 MeV. (27) Our results for the N f = 2 massless case are shown in Fig. 8, in which we show both the quark condensate and the eigenvalue gap values extrapolated to the infinite volume limit. Both of them indicate the same location of the critical temperature, which is finally determined by a fit with expected critical exponent (solid curve). Translating our scale to approximate absolute temperature, we found that our simulations cover the temperature range of about 80 MeV < T < 150 MeV.  (26) and is a fit to the four data points just below Tc. The overall normalization of Σ(T ) is arbitrary.

C. Effects of nonzero quark mass
A nonzero quark mass explicitly breaks chiral symmetry in the QCD Lagrangian. As seen in the previouslymentioned lattice works, the dynamics of the chiral phase transition in particular are sensitive to the masses of the lightest quark flavors. A full treatment of a theory with nonzero quark masses would require the following modifications to our dyon partition function: (i) A generalization of the perturbative quark potential (8) to arbitrary quark mass, the form of which is given in Ref. [24]. (ii) A generalization of the hopping matrix elements T ij to arbitrary quark mass, which is not yet known. (iii) The inclusion of quark mass terms on the diagonal elements of the hopping matrix as shown in Eq. (17). We do not do such a full treatment in this work. Instead we include only the quark mass term on the diagonals of the hopping matrix to demonstrate the qualitative impact it has on the eigenvalue spectrum. The nonzero quark mass effectively adds new diagrams to the fermionic determinant in which single dyons are allowed to have closed loops with a mass insertion that flips the chirality of the quark. The mass mediates the behavior of the quark-induced potential, driving it to a constant value at large distances, rather than remaining linear as in the massless case.
Eigenvalue spectra for a small, nonzero quark mass are compared to the massless case in Fig. 9. The nonzero quark mass allows for near-zero eigenvalues even at finite volume. The quark mass smooths the distribution of eigenvalues in the vicinity λ ∼ m. At eigenvalues λ > m, the distributions are the same as in the massless case. In the broken phase, increasing the quark mass reduces the value of the condensate. It is known for example that the strange quark condensate is smaller than the up quark condensate [43]. The increased mass causes the condensate to decrease more slowly, becoming an analytic crossover with increased (now psuedocritical) T c . One can see from Fig. 9 (right) that there is a nonzero condensate in the nonzero mass case, while it has already reduced to zero in the massless case. The finite-mass condensate never exactly goes to zero, as there is no longer an exact symmetry.

VII. SUMMARY AND DISCUSSION
In this work we have performed numerical simulations of a semiclassical ensemble of SU (3) instanton dyons with N f = 2 flavors of massless quarks. Integration over the dyon degrees of freedom was performed in a periodic 3D box via Monte-Carlo methods. From tens of thousands of simulations, the properties of the ensemble in the thermodynamic limit were determined.
The main addition to the simulations stemming from the inclusion of quarks is the computation of the fermionic determinant. This is approximated by the socalled hopping matrix, which contains only the subspace spanned by the quark zero modes. Quarks then 'hop' between L andL dyons generating an effective interaction between all such dyons in the ensemble. These interactions, which can be dominated by single LL pairs or collective modes involving the overlap of many zero modes (see Fig. 2), determine the state of the chiral symmetry breaking.
The two phase transitions -deconfinement and chiral symmetry restoration -were observed. Confinement is studied by the value of the average Polyakov loop P (T ) , seen in Fig. 4. Indeed we find that the inclusion of quarks changes the deconfinement transition from first order to one compatible with that of the second-order transition in the O(4) universality class.
The near-zero-mode zone of the Dirac eigenvalue spectrum (Fig. 6) is used to determine the zero-eigenvalue density. This is directly related to the quark condensate via the Banks-Casher relation (23). We performed simulations at three different ensemble sizes in order to observe the non-trivial volume dependence of the spec-tra (Fig. 7) and measure both the quark condensate and the eigenvalue gap as functions of the temperature. Fig.  8 shows that both observables see nearly-simultaneous transitions to/from zero giving consistent determinations of the critical temperature T c . Finally, we show that a small, but nonzero quark mass smooths the distribution and produces near-zero eigenvalues at higher temperatures, increasing T c .
Our results suggest that both phase transitions are driven primarily by the dyon densities. Confinement requires a sufficient density of dyons to overcome the perturbative quark and gluon interactions and shift the minimum to the confining holonomy. Chiral symmetry is similarly broken by a large density of L andL dyons causing significant overlap between zero modes, leading to a dominance of large quark hopping loops, producing near-zero eigenvalues.
Let us conclude with some discussion of the dyon model itself. It should be reminded that, as with the pure SU (3) theory, our model contains parameters related to the short-range classical interactions (namely V 0 and x 0 ) which are phenomenological choices not known from first principles. With the hopping matrix elements we choose a simple parameterization that does not include the effects of interference from nearby M i dyons. Both of these aspects of the model should be studied more rigorously in order to improve the quantitative predictions of the model.
While we have been able to identify both phase transition, compared with the pure SU (3) work [10] the lack of a jump in the order parameters makes precise determinations of the critical temperatures more difficult. The deconfinement transition is slow and fits to potential forms yield a large variance in T deconf . With the chiral transition, the drop in the condensate is much sharper, but requires better control of finite volume effects and different methods of extrapolation can modestly modify the determination of T c (see Appendix B).
Despite its relative simplicity, there are some advantages to the dyon model compared with analogous lattice studies. The most obvious is the computational cost. Our largest simulations with N D 360 involve integration over ∼ 1000 degrees of freedom, while modern lattice simulations can involve some hundreds of millions of degrees of freedom. Additionally our largest simulations contain O(100) instantons at a time, significantly more than the number of instantons that are observed on a single time slice of lattice simulations. Thus, our simulations correspond to much larger spatial volumes than are used on the lattice. Lastly we are able to work directly with massless quarks, where the recent lattice results [28][29][30][31][32][33] require nonzero quark masses and an extrapolation to the chiral limit.
The quark zero mode is localized on the dyon i such that µ i < z < µ i+1 , which is the L dyon associated with i = 3.
Taking the necessary derivatives to write out either density directly in terms of the dyons' coordinates, even after simplification, results in formulae which are dozens of lines long as Mathematica outputs. These are far too complicated to compute at each step of a Metropolis update, hence the use of a simple parameterization in this work. In the limit that the dyons are well separated, one can see that the general zero-mode density given here is exactly what is predicted by Eq. (19), up to an overall normalization constant.

Appendix B: Comparison of Infinite-Volume Extrapolations
The determination of the chiral condensate depends a choice of interpolation function between different system sizes (25). Rather than interpolating between all three volumes, one could consider using just two of the volumes. In particular we consider scaling directly between volumes V and 3V and scaling between 2V and 3V with the functions As with the interpolating function used in the main text, we set a maximum value of 1 on the scaling factors. Each of these functions results in different results for the chiral condensate and are compared in Fig. 10. There are a few qualitative differences between them. For several of the data points, all three give the exact same answers as the range of the suppressed region scales faster than 1/V and the scaling functions are all taking on the ceiling value of 1. At the lowest values the function using all three volumes decreases slowly rather than increasing like the other two. This is because the scaling factor from V to 2V is less than 1 while the others are 1. This may be an indication that at these temperatures, where the density is highest and the volumes are smallest, the system may not be large enough for the mesoscopic scaling relations (24) to apply. One should note that it is not the total number of dyons N D , but the number of zero modes 2N L which is relevant to the eigenvalue distributions. For the smallest ensemble sizes, 2N L ∼ 20 − 25.
From Fig. 10 one can see that the value of T c can vary by about 10%, depending on the choice of interpolating function. Of course, one could also consider other functions which may produce even more varied results.
One could also consider using only two ensemble sizes for the linear fit to the eigenvalue gap data. Doing so however, one finds much less dependence on the choice of which sizes to include. At low T the gaps also show better agreement with the expected 1/V scaling and do not need an interpolating function. Thus we conclude that the eigenvalue gap ∆(T ) is a more stable and reliable indicator of which chiral phase the system is in. Of course, the best way to improve the results for either observable is to continue to go to larger volumes.