  Topological phenomena in gauge theories have long been recognized as the
driving force for chiral symmetry breaking and confinement. These phenomena can
be conveniently investigated in the semi-classical picture, in which the
topological charge is entirely carried by (anti-)self-dual gauge
configurations. In such an approach, it has been shown that near the critical
temperature, the non-zero expectation value of the Polyakov loop (holonomy)
triggers the "Higgsing" of the color group, generating the splitting of
instantons into $N_c$ self-dual dyons.
  A number of lattice simulations have provided some evidence for such dyons,
and traced their relation with specific observables, such as the Dirac
eigenvalue spectrum. In this work, we formulate a model, based on one-loop
partition function and including Coulomb interaction, screening and fermion zee
modes. We then perform the first numerical Monte Carlo simulations of a
statistical ensemble of self-dual dyons,as a function of their density, quark
mass and the number of flavors. We study different dyonic two-point correlation
functions and we compute the Dirac spectrum, as a function of the ensemble
diluteness and the number of quark flavors.
