1-\cos\left(\frac{2\pi n_{\mu\nu}}{N}\right)\right\}
   \simeq 2\pi^2\beta\sum_{\mu\neq\nu}n_{\mu\nu}^2. 
\label{EQ:potential_difference}\end{aligned}$$ For the generic twist, it becomes $$\Delta S=
\begin{cases}
24\pi^2\beta m^2L^2\qquad & \mbox{(symmetric form)},\\
8\pi^2\beta m^2L^2\qquad  & \mbox{(skew-diagonal form)}.
\label{EQ:potential_difference_generic}
\end{cases}$$ Note that the symmetric form is roughly three times more stable than the skew-diagonal form if both twists have equal quantum fluctuations. This is the reason that the $\mathbb{Z}_N^4$ symmetry breaking for the skew-diagonal form can occur at quite smaller $N$ than that for the symmetric form, as is observed in our simulation.

### Quantum fluctuations and symmetry breaking {#SSEC:fluctuation .unnumbered}

Going away from the weak coupling limit, the system has quantum fluctuations. Here we naively expect that the $\mathbb{Z}_N^4$ symmetry is broken if the fluctuation around twist-eater configuration exceeds the energy difference $\Delta S$. Because the system describes $O(N^2)$ interacting gluons, it is natural to assume that their quantum fluctuations provide $O(N^2)$ value to the effective action. For the generic twist the quantum fluctuation is $O(m^2L^4)$. Combined with the fact (\[EQ:potential\_difference\_generic\]), we can estimate the critical point $\beta_c^L$ as $$\beta_c^L\sim L^2,\label{EQ:linear_behavior}$$ which is consistent with the numerical results (\[EQ:bc\_L\_minimal\_sym\]), (\[EQ:bc\_L\_minimal\_skew\]) and (\[EQ:bc\_L\_generic\_skew\]). In addition we can explain the difference of the coefficient of $L^2$ in (\[EQ:bc\_L\_minimal\_sym\]), (\[EQ:bc\_L\_minimal\_skew\]) and (\[EQ:bc\_L\_generic\_skew\]) between the symmetric and 