  We study $Z_N$ symmetry in $SU(N)$ gauge theories in the presence of matter
fields in the fundamental representation. The Polyakov loop dependence of the
free energy is calculated by restricting the lattice partition function
integration to matter fields which are uniform in spatial directions and gauge
fields with vanishing spatial components. We show that in the limit of large
number of temporal sites the explicit breaking of $Z_N$ symmetry vanishes,
driven by dominance of the density of states. We argue that the spatial links
as well as the spatial modes of the matter fields determine the boundaries
separating regions where $Z_N$ symmetry is realised from the rest.
