ambiguous because the breakdown of the $\mathbb{Z}_N^4$ symmetry is a first-order transition. As an order parameter for detecting the $\mathbb{Z}_N^4$ breakdown, we measure the Polyakov lines $$\begin{aligned}
P_\mu\equiv\left\langle\left|\frac{1}{N}{\rm Tr}~U_{\mu}\right|\right\rangle.\end{aligned}$$

Simulation results
------------------

### Minimal symmetric twist {#minimal-symmetric-twist .unnumbered}

First of all we treat the minimal symmetric twist (\[EQ:minimal\_sym\_form\]). This twist is the most standard and is also used in the paper [@Teper:2006sp]. In our study we only investigate the first $\mathbb{Z}_N^4$ symmetry breaking point from weak coupling limit, that is, $\beta_c^L$ for this twist. (For more detailed information about the symmetry breaking phenomena, see [@Teper:2006sp].) The obtained results are in table \[TAB:minimal symmetric twist\] and plotted in figure \[FIG:minimal symmetric twist\]. The symmetry breaking points and patterns ($\mathbb{Z}_N^4\xrightarrow{\beta_c^L}\mathbb{Z}_N^3$ for $N=100$; $\mathbb{Z}_N^4\xrightarrow{\beta_c^L}\mathbb{Z}_N^0$ for $N>100$) are consistent with the results in [@Teper:2006sp] up to $N=144$. In this work we explore the simulation for larger $N$. From figure \[FIG:minimal symmetric twist\] we can find clear linear dependence of $\beta_c^L$ on $N(=L^2)$ for $N\gtrsim169$. The fitted result in linear function using $N\geq169$ data is $$\beta_c^L\sim 0.0011N+0.21.
\label{EQ:bc_L_minimal_sym}$$ A theoretical argument for this linear behavior is discussed in section \[SEC:discussions\].

  $N$     $L$        $\beta_c^L$      $N$     $L$        $\beta_c^L$
  ------- ------ -------------------- ------- ------ --------------------
  $100$   $10$    $0.3525\pm 