Lattice ${\mathbb C} P^{N-1}$ model with ${\mathbb Z}_{N}$ twisted boundary condition: bions, adiabatic continuity and pseudo-entropy

We investigate the lattice ${\mathbb C} P^{N-1}$ sigma model on $S_{s}^{1}$(large) $\times$ $S_{\tau}^{1}$(small) with the ${\mathbb Z}_{N}$ symmetric twisted boundary condition, where a sufficiently large ratio of the circumferences ($L_{s}\gg L_{\tau}$) is taken to approximate ${\mathbb R} \times S^1$. We find that the expectation value of the Polyakov loop, which is an order parameter of the ${\mathbb Z}_N$ symmetry, remains consistent with zero ($|\langle P\rangle|\sim 0$) from small to relatively large inverse coupling $\beta$ (from large to small $L_{\tau}$). As $\beta$ increases, the distribution of the Polyakov loop on the complex plane, which concentrates around the origin for small $\beta$, isotropically spreads and forms a regular $N$-sided-polygon shape (e.g. pentagon for $N=5$), leading to $|\langle P\rangle| \sim 0$. By investigating the dependence of the Polyakov loop on $S_{s}^{1}$ direction, we also verify the existence of fractional instantons and bions, which cause tunneling transition between the classical $N$ vacua and stabilize the ${\mathbb Z}_{N}$ symmetry. Even for quite high $\beta$, we find that a regular-polygon shape of the Polyakov-loop distribution, even if it is broken, tends to be restored and $|\langle P\rangle|$ gets smaller as the number of samples increases. To discuss the adiabatic continuity of the vacuum structure from another viewpoint, we calculate the $\beta$ dependence of ``pseudo-entropy"density $\propto\langle T_{xx}-T_{\tau\tau}\rangle$. The result is consistent with the absence of a phase transition between large and small $\beta$ regions.

The model on R × S 1 with Z N symmetric twisted boundary conditions (Z N -TBC) also attracts a lot of attention since the Z N symmetry, whose order parameter is the expectation value of the Polyakov loop operator, is exact in this model. This model admits fractional instantons [10,[50][51][52][53], namely instantons with fractionally quantized topological charge, typically 1/N quantized one for the Z N -TBC (see Refs. [54][55][56] for fractional instantons in other models with TBC). Then, it has been conjectured [57][58][59][60] that the Z N -symmetric vacuum of the model on R × S 1 is continuously connected to that on R 2 due to the tunneling transition by fractional instantons. If there exists such an adiabatic continuity of the vacuum, it gives us a deeper understanding of the model on R 2 through the weak coupling analysis on R × S 1 . In particular, the nontrivial relation between perturbative and nonperturbative contributions, called "the resurgent structure" [61][62][63][64][65][66][67][68][69][70][71], which has been intensively studied in this model [57,58,[72][73][74][75][76][77][78][79][80][81][82][83][84][85], is expected to play an important role in understanding the relation between the weak and strong coupling regime of the model. From this perspective, it is of great importance to study the CP N −1 model on R × S 1 with the Z N symmetric twisted boundary condition, with particular attentions to the Z N symmetry and its order parameter, namely the expectation value of the Polyakov loop P .
In this paper, we investigate the CP N −1 model on S 1 s (large) × S 1 τ (small) with the Z N symmetric twisted boundary condition by lattice Monte Carlo simulations. The ratio of the circumferences L s /L τ is taken to be sufficiently large so that the model on S 1 s ×S 1 τ approximately describes that on R × S 1 τ . We focus on the distributions and expectation values of the Polyakov loop, the dependence of the argument of the Polyakov loop on L s . In addition, we also investigate the "pseudo-entropy" density ∝ T xx −T τ τ , which is the counterpart of the thermal entropy density in the Z N -TBC case.
For small inverse coupling β (or large L τ ), the value of the Polyakov loop for each configuration concentrates around the origin on the complex plane and its expectation value is zero (| P | = 0) as a result of the exact Z N symmetry. The distribution of the Polyakov loop spread around the origin broadly and isotropically as β increases. We find that, in high β regions (or small L τ regions), the distribution forms a regular N -sided-polygon shape and the expectation value of the Polyakov loop is still consistent with zero (| P | ∼ 0). By studying the L s -dependence of the Polyakov loop, we also show the existence of fractional instantons and bions, which cause tunneling transition among the classical N vacua, leading to the stabilization of the Z N -symmetric vacuum. For much higher β, the statistics in our simulation are less than the auto-correlation time for | P | and the results of the simulation gets unreliable. Even in such high β regions, we find that a regular N -sided-polygon shape of the Polyakov-loop distribution tends to be restored and | P | gets smaller by increasing the number of samples. These results indicate the stability of the Z N symmetry and also imply that the seemingly broken Z N symmetry at extremely high β can be an artifact due to insufficient statistics. In order to study the adiabatic continuity from another viewpoint, we calculate the β dependence of the pseudo-entropy density, s ∝ T xx − T τ τ , and find no phase transition as β increases from small to large values. Furthermore, our result indicates that the pseudo-entropy density vanishes in the large-N limit. It suggests the volume independence in the whole β regime in the large-N limit, consistent with the analytical study of "the large-N volume independence" in Ref. [60]. This paper is organized as follows. In Sec. II, we introduce the model in the continuum limit and review its properties. In Sec III, we review the results of lattice simulation for the model with the PBC. In Sec. IV, we show the results of the lattice Monte Carlo simulation for the model with the Z N TBC. In Sec. V, we measure the pseudo-entropy density and discuss its implication on the adiabatic continuity. Section VI is devoted to a summary and discussion.
II. TWO-DIMENSIONAL CP N −1 MODEL In this paper, we investigate CP N −1 = SU(N )/(SU(N − 1) × U(1)) sigma models (without the topological θ-term) on R × S 1 , with attentions to the Z N -TBC, 1/N fractional instantons and Z N intertwined symmetry. In this section we review these notions.

A. Basics of the model
Let ω(x) be an N -component vector of complex scalar fields, and φ(x) be a normalized complex The action of the CP N −1 model in Euclidean two dimensions is given by where g 0 is the bare coupling constant, d 2 x ≡ dxdτ and the indices µ, ν = 1, 2 label the x, τ directions. The covariant derivative is defined as D µ = ∂ µ + iA µ with a composite gauge field where the Z N center is removed since it coincides with a subgroup of the U(1) gauge symmetry and hence is redundant. All through this paper, we consider (or approximate) the model on R × S 1 , and regard x and τ as coordinates of uncompactified (large) and compactified (small) directions, respectively.
This model has instanton solutions characterized by the topological charge representing The simplest case, or the CP 1 model, is equivalent to the O(3) nonlinear sigma model, thus it can be also described by three real scalar fields m(x) = (m 1 (x), m 2 (x), m 3 (x)) with a constraint m(x) 2 = 1. Its action is given by where the relation between the real scalar field m(x) and the complex two-component complex with the Pauli matrices σ. In this description, the configuration φ = (1, 0) T corresponds to the north pole m = (0, 0, 1), while φ = (0, 1) T corresponds to the south pole m = (0, 0, −1).
The boundary condition we mainly focus on is the Z N -symmetric twisted boundary condition (Z N TBC). The Z N -TBC in a compactified direction can be expressed as where L τ is the compactification circumference. We note that this Z N -TBC is equivalent to the existence of the following Wilson-loop holonomy of the background SU(N ) gauge field in the compactified direction: Let us discuss the symmetry of the present model. Consider the Z N subgroup of the flavor SU(N ) transformation generated by This transformation alone does not keep the Z N -TBC but becomes a symmetry which does not change the boundary condition if it is combined with the large gauge transformation Explicitly, we can confirm that the intertwined transformation does not change the boundary condition as where we have used the fact that the matrix Ω is invariant under the combined transformation Therefore thus model has the intertwined Z N symmetry, which is the invariance under the simul- We now introduce fractional instantons [10,50,51] in the CP 1 model satisfying the Z 2 -TBC, which is given by or equivalently By defining the complex coordinate z = x + iτ , we can express one of the fractional instanton solutions as where two real moduli parameters (a position and phase) are combined into a complex constant a. x 0 = 1 π log 1 |a| . Under the Z 2 -TBC, the U (1) modulus is twisted by Z 2 along the domain wall. It means that, when the constant τ slice is changed from τ = 0 to τ = L τ , the path in the target space sweeps a half of the sphere of the target space. Therefore, the configuration is a map from the space R × S 1 to a half of the target space. This is the reason why it carries a half of the unit instanton charge Q = 1/2 while the anti-BPS configurations carry Q = −1/2. One finds that the topological charge density has no dependence on the compactified direction τ .
The fractional instantons in the CP N −1 model are classified in a parallel manner. The configuration (14) of the CP 1 model can be generalized to the N -vector ω for the CP N −1 model with the These configurations carry 1/N unit of the instanton charge (Q = 1/N ). These fractional instantons correspond to the tunneling transitions (domain walls) among N classical vacua. It means that they stabilize the Z N -symmetric vacua at the quantum level. In the quantum-mechanical limit with L τ → 0, we can easily show that the Z N symmetry is preserved at the quantum level due to the fractional instantons. The conjecture on the adiabatic continuity of the vacuum structure states [60] that "the Z N -symmetric vacua continue to be unbroken from the small S 1 regime to the large S 1 regime". One of the purposes of the present work is to study the conjecture in terms of lattice Monte Carlo simulation.

A. Lattice setup
We now review the results for PBC, in particular the expectation value of the Polyakov loop [47]. The lattice action of the two-dimensional CP N −1 sigma model [36-39, 41, 44, 47] is withφ n · φ n = 1 and λ n,µ being a link variable corresponding to the auxiliary U(1) gauge field.
The lattice sites n = (n x , n τ ) run as n x = 1, · · · , N s and n τ = 1, · · · , N τ and N β is equivalent to the inverse of the bare coupling 1/g 2 0 . The circumferences of S 1 s and S 1 τ are given in terms of the lattice spacing a as L s = N s a and L τ = N τ a, respectively.
The renormalization group tells us the relation between β and the lattice spacing a where the renormalized coupling in the M S scheme diverges at Λ M S . By comparing Λ M S [37] and the lattice Λ scale Λ lat for Eq. (16), we obtain the following relation, which relates a with β for each N , where Λ lat works as a reference scale. The parameter set with L s L τ was taken to approximately simulate the model on R × S 1 , where L τ corresponds to an inverse temperature 1/T and smaller L τ or higher β with fixed N τ corresponds to higher T .

B. Lattice simulation for periodic boundary condition
The Polyakov loop on the lattice is given by It is an order parameter of the Z N symmetry, whose transformation is given by While this is an exact symmetry and the expectation value of the Polyakov loop P is its order parameter in the model with Z N -TBC, the model with the PBC is not Z N -symmetric. However, P is still of physical importance in the PBC case. Based on the relation between the Polyakov and Wilson loops via the clustering property, one finds that the expectation value of the Polyakov loop vanishes ( P = 0) in the confinement phase with a nonzero string tension in the L τ L s system.
Thus, P can be adopted as the order parameter of the confinement-deconfinement transition in the system. Below, we summarize the main results of the lattice simulation for the model with PBC. All these results strongly indicate a crossover behavior between the confinement and deconfinement phases.
1. The distribution plots of the Polyakov loop for N = 3, 5, 10, 20 exhibit that the Polyakov loops are distributed around the origin at low β and gradually get away from the origin at higher β. of | P | is independent of β in β β c since the IR scale of the system is the confinement scale, but not 1/L τ , in the confining phase. We here denote β c as a pseudo-critical β, which is defined from the peak position of χ |P | = V ( |P | 2 − |P | 2 ). As β increases across β c , | P | gradually increases and we find | P | = 0 for high β (small L τ ).
3. The Polyakov-loop susceptibility χ |P | as a function of 1/L τ shows that the peak is broad for small N but gets sharper as N increases.
4. The volume dependence of the peak value of χ |P | is measured by simulating N s = 40, 80, 120, 160, 200 with N τ = 8 fixed. The fit of data points for N s = 80, 120, 160, 200 by a function χ |P | ,max = a + cN p s [86] indicates that the best fit values of the exponent are p < 1 for N = 3, 5, 10, 20, in conformity with a crossover behavior.
This Z N symmetry is conjectured to be unbroken even at high temperature (small L τ ) due to transitions among the Z N vacua via the fractional instantons [57,58].
We here perform the Monte Carlo simulation with the Z N -TBC imposed on the small com- The former, | P |, is a genuine order parameter of the Z N symmetry, while the latter, |P | , is not an order parameter of Z N symmetry, indicating the averaged distance of the Polyakov-loop values from the origin in the complex plane. For the case of the PBC, these two quantities show a crossover behavior and have an identical behavior except for a small β region, where we have | P | = 0 while |P | takes a finite value due to the finite volume effect [47]. We show below that it is not the case with the Z N -TBC. the bin size N bin to be larger than the auto-correlation time. We also show the distribution plots of the Polyakov loop for these cases with several different β in Fig. 2.
In Fig. 1, |P | (blue points) gradually increases in β c < β for each N , e.g. β c ≈ 1.0 for N = 3 and β c ≈ 0.5 for N = 20, since the distribution of the Polyakov loop gradually spreads as shown in Fig. 2. We here define β c as the value of β at which |P | starts to increase and 1/L τ becomes the IR scale of the system. In β < β c regime, the values of |P | is independent of β and take the same values as those for PBC as we will show in Fig. 8. It indicates that the IR scale of the system changes from Λ CP N −1 to 1/L τ at β c .
On the other hand, the genuine Z N order parameter | P | (red points) is almost consistent to zero within statistical errors even at a high β region (β > β c ). This behavior originates in the fact that the distribution of the Polyakov loop spreads isotropically around the origin. Let us investigate this Z N order parameter | P | in more detail by comparing it to the distribution plot in Fig. 2 We now discuss the validity of our simulation results. As will be discussed in App. A, we have estimated the auto-correlation times for | P | by studying the bin-size dependence of the Jackknife error bars. Based on the estimated auto-correlation time, we consider that our simulations, whose main results are shown in Figs. 1 and 2, are performed with the larger number of samples than the auto-correlation times. On the other hand, for much higher β, the auto-correlation times seem to be larger than the number of samples (800, 000 for N = 3, 5 and 600, 000 for N = 10, 20), and thus we cannot obtain meaningful results of the Polyakov-loop expectation values in the β region so far. Indeed, in this β region, the regular N -sided polygon shapes of the Polyakov-loop distribution lose their shapes as shown in Fig. 3. We consider that this behavior is just an artifact due to the shortage of statistics since the regular N -sided polygon shapes of the Polyakov-loop distribution are getting restored as we adopt a larger number of samples even for these cases as we will discuss in the next-next subsection.

B. Fractional instantons
We next focus on each of the configurations constituting the N -sided polygon-shaped distributions in Fig. 2. As we have discussed in Sec. II, the fractional instantons work to stabilize the n x dependence of arg[P (n x )], which is defined as n x is the lattice coordinate for S 1 s (large) direction. We note that arg[P (n x )] ≈ A τ dτ describes the vacuum transition process since the topological charge Q for R × S 1 is given by [51] For our lattice setup on S 1 s (large) × S 1 τ (small), the PBC is also imposed on n x direction, thus the total Q should be an integer [89]. However, the transition process can be nontrivial, and it can include nontrivial configurations such as bions.
We now pick up two of configurations (A) and (B) constituting the polygon-shaped distribution at β = 1.6 for N = 3 in Fig. 4(Left). In Fig. 4(Center), we depict arg[P (n x )] for the configura- Since the relatively low-β configurations suffer from large fluctuations, it is not easy to identify fractional instanton configurations without the auxiliary lines as shown in Fig. 4

(Center)(Right).
For high-β, these configurations, if exist, become clean because of less fluctuations. Although the polygon-shaped distribution tends to be broken at extremely high β due to the shortage of statistics as we have discussed, the fractional instanton configuration can be found at certain probability even at quite high β. For instance, we find the fractional instanton configurations at ∼ 10% probability at β = 3.0 for N = 3. We depict arg[P (n x )] for two of such configurations at β = 3.0 in Fig. 5    We now make comments on the reason why the auto-correlation time increases and why we need larger statistics for higher β. The theoretical reason is a rapid decrease of the transition probability among the N classical vacua. The transition probability between adjacent N classical vacua due to 1/N fractional instantons with the action S = 2π N g 2 is expressed as Although g 2 should be the renormalized coupling, we can approximate it by the bare coupling g 2 0 = 1 N β in the weak-coupling limit (corresponding to small L τ or large β limit). Therefore, the transition probability exponentially decreases when β increases. In the Monte Carlo simulation, it means that we need exponentially large N sweep to observe a true quantum vacuum, which is expected to have a small expectation value of the Polyakov loop (| P | ∼ 0). This is one of possible reasons why the auto-correlation time surges for larger β. It is also notable that the vacuum "freezes" for very high β and the auto-correlation time seems to get small in appearance because of the small transition rate among the classical vacua. However it is as an artifact and a large autocorrelation time should be observed if one adopts exponentially large statistics. In other words, the statistical errors of | P | in extremely high β regions may be underestimated in the simulation due to the freezing of the vacuum.
Since we do not have enough statistics at β > 2.0 for any N in our simulations, we cannot have a strong conclusion on the adiabatic continuity of the Z N -symmetry for the whole β regime. However, we emphasize that the Z N -symmetric regions exist even in β c < β, where the IR scale is no longer an ordinary confinement scale as shown in Sec. IV A. Furthermore, the results in this subsection imply that we may be able to observe the adiabatic continuity by adopting an exponentially large number of samples N sweep . This speculation will be also supported by the arguments in Sec. IV E. transitions (or the absence of sudden expanses of the distribution) in the Z N -twisted model.
In the previous work [47], it was shown that the β dependence of |P | undergoes a crossover behavior in the CP N −1 model with PBC. We here make use of this fact and make comparison between |P | of the Z N -TBC and PBC. If we find gentler increase of |P | in the Z N -TBC than that in the PBC as β increases, we can conclude the absence of phase transitions between large-L τ and small-L τ regions at least for this quantity.
We depict |P | for the Z N -TBC and PBC in Fig. 8. The reason why we plot them from small to quite large β (β < 3.0) is that the auto-correlation time for this quantity is much smaller than that for | P | and the results of the simulations are reliable. It is obvious that, as β increases, |P | for the Z N -TBC gets larger more gradually than that for the PBC. It means that |P | for the Z N -TBC undergoes a crossover behavior between large-L τ and small-L τ regions, thus there are no phase transitions regarding |P | in this model. We consider that this is an affirmative fact the absence of sudden change of vacuum structure.

E. Large-volume simulation
In this subsection, we show results of the simulation with larger volume (N s , N τ ) = (400, 12) for N = 3. In particular, we perform an exploratory simulation at very high β as β = 4.0 with N sweep = 198, 000 and obtain distribution plots of the Polyakov loop. Since an auto-correlation time for | P | at β = 4.0 is speculated to be quite large, this N sweep is probably insufficient. Therefore the reliability of results in this simulation is lower than those in other parts of this paper, and we do not show a result of | P | but just discuss one of distribution plots. This is why we refer to the simulation in this subsection as an "exploratory" one.
In Fig. 9(Left), we depict one of the distributions plots of the Polyakov loop, which has a partially broken regular-triangle shape. We check that the history of arg[P ] for this distribution is quite random in Fig. 9(Right), implying that the transition among the Z N vacua occurs perpetually.
This result implies that, even at quite high β, the regular-polygon distribution of the Polyakov loop appears with a certain probability (∼ 5%) at least for a larger volume (N s , N τ ) = (400, 12).
In this simulation, the fractional instanton configurations are also observed. We pick up two con- entropy density" for Z N -TBC, since it still carries similar properties to the thermal entropy. In terms of the energy-momentum tensor (EMT) T xx , T τ τ , the pseudo-entropy density s in the large volume limit is given by where the Z N -TBC along the τ direction is imposed. On the lattice, EMT is defined [90] as where the vacuum expectation value of the trace part is subtracted [91,92]. We will adopt the bare coupling constant to calculate these quantities since it well approximates the renormalized coupling in the weak coupling regime.
Let us discuss the simulation results of the pseudo-entropy density s/(N T ). The β dependences of the pseudo-entropy density for N = 3, 5, 10, 20 are depicted in Fig. 11(left). The pseudo-entropy density becomes non-zero around a certain β and monotonically decreases. In the high-β regime, the which is depicted as the broken curve. For small N such as N = 3, 5 in the figure, small deviations from the analytical result are found. We speculate that it indicates 1/N 2 corrections to the leadingorder results.
It is worth noting that, in the large-N limit, the pseudo-entropy in the high β regime seems to go to zero. Thus, the difference between T xx and T τ τ disappears even though the anisotropy between x and τ directions is large. It suggests the volume independence of the EMT and the action density.
This speculated property in the large-N limit is consistent with the analytical result proposed in It is also notable that the pseudo-entropy densities do not exhibit any special behaviors such as phase transition at high β. We consider that these results of the pseudo-entropy density also support the adiabatic continuity of the model between R × S 1 and R 2 although this quantity is not sensitive to the breaking of the Z N symmetry.
In the end of this section, we make a brief comment on the pseudo-entropy in the QCD-like models. The pseudo-entropy can be also calculated in the lattice simulation of N -flavor QCD with the Z N twisted boundary condition (Z N -QCD) [87,88]. The question whether it is negative or not in Z N -QCD is quite nontrivial and will be addressed in future works.

VI. SUMMARY AND DISCUSSION
In this paper, we have reported the Monte Carlo simulations for the CP N −1 model on S 1 s (large)× S 1 τ (small) with the Z N -TBC, where sufficiently large ratio of the circumferences is taken to approximate the model on R × S 1 . We have found that the β dependence of expectation values of the Polyakov loop differs significantly from those for the PBC. Our results so far seem to give evidences in support for the conjecture of the adiabatic continuity of the vacuum structure. We here summarize our main results: 1. By studying |P | , we find that the Z N -twisted model has the characteristic β denoted as β c , at which the IR scale of the system changes from Λ CP N −1 to 1/L τ .
2. The order parameter of the Z N symmetry, | P | ∼ 0, continues to be consistent with zero in the both β < β c and β > β c regions. It means that the Z N symmetry is unbroken even in the high β regime. The distribution of the Polyakov loop in β > β c forms regular N -sided-polygon shapes.
3. By investigating the dependence of the Polyakov-loop phase on the coordinate of the S 1 s direction, the existence of fractional instantons and bions causing the transition among the classical N vacua is verified. We argue that fractional instantons work to stabilize the Z N symmetry.
4. The result on | P | is unreliable for quite high β since the statistics we adopt in the simulation are less than the auto-correlation time. Even for such a high β, a regular N -sided-polygon shape of the Polyakov-loop distribution tends to be restored and | P | gets smaller by increasing the number of samples. 5. In a larger-volume simulation, the N -sided polygon-shape in the Polyakov loop distribution appears even for extremely high β with a certain probability. We note that it is an exploratory calculation and the number of samples is not sufficient. 6. The β dependence of the pseudo-entropy density s = T xx − T τ τ /T implies the absence of phase transitions between the large and small circumference regions. Furthermore, the pseudo-entropy density vanishes in the large-N limit. It is consistent with the volume independence in the whole β regime in the large-N limit.
This work helps deepen the understanding on the symmetry and phase diagram of the CP N −1 model on the compactified spacetime with the Z N -TBC. Our main focus was the adiabatic continuity of the Z N symmetric phase, which is essential to apply the resurgence theory to the model on R × S 1 in the decompactified limit. Although our results are not conclusive to the conjecture of the adiabatic continuity in the model, there are various implications for future studies on the topic.
In the end of this paper, let us introduce our next plan to judge the adiabatic continuity of the Z N symmetry in the model: (1) At quite high β (small L τ ), we generate configurations which give polygon-shape distribution and a small expectation value of Polyakov loop | P | ∼ 0 (e.g. the case of Fig. 9).
(2) We pick up one of the above configurations and use it as an initial configuration to generate configurations at slightly lower β.
(3) We repeat this procedure and generate configurations from high to low β.
(4) We investigate whether | P | ∼ 0 continues form high to low β (from small to large L τ ).
The point is that the generation of configurations in this proposal starts from high β while we started it from low β in the simulation of the present work. We can judge the validity of the adiabatic continuity through this procedure in principle, although we need to realize the "adiabatic" decrease of β which is quite hard to achieve in lattice simulations. In Fig. 12, we depict it for N = 3 and β = 1.8, 1.9. For β = 1.8, the plateau appears around N bin = 160000. For β = 1.9, the plateau does not appear within the number of samples we adopt in this work. Since the value of N bin at which the plateau appears corresponds to auto-correlation time, we roughly estimate that the simulations for β ≤ 1.8 are performed with sufficient statistics beyond the auto-correlation times. For higher β, we do not reach a plateau within the statistics, thus we consider that they are less than the auto-correlation time and the simulation results are not reliable.
In Fig. 13, we depict the bin-size dependence of errors of the Polyakov-loop expectation values for N = 5 and β = 1.6, 1.8. For β = 1.6, the plateau appears around N bin = 160000. For β = 1.8, the plateau does not appear within the number of samples we adopt in this work. (For β = 1.7, we cannot judge whether or not the plateau exists due to the fluctuation.) As with the case of N = 3, we roughly estimate the simulations for β ≤ 1.6 are performed with sufficient number of samples beyond the auto-correlation times while they are insufficient for higher β.
We have done the same analysis for N = 10, 20. We roughly estimate that the simulations for  Appendix B: Thermal entropy density from the free energy density In this appendix, we derive thermal entropy density for a single free scalar field with the TBC.
The TBC for a scalar field φ is given by Then, the TBC analog of the partition function is given by the the TBC analog of the free energy F = − 1 Lτ log Z can be rewritten as F = L s d 2 p (2π) 2 e i(nτ (Lτ pτ −α)+nsLsps) log(p 2 + m 2 ) + const , where denotes the summation over (n τ , n s ) ∈ Z 2 excluding the term with (n τ , n s ) = (0, 0) (denoted as const). Performing the integration, we obtain F = − e −inτ α mL s π L 2 τ n 2 τ + L 2 s n 2 s K 1 (m L 2 τ n 2 τ + L 2 s n 2 s ) , where K 1 stands for the modified Bessel function of the second kind. In the large L s regime, the leading contribution of the free energy density is given by the terms with n s = 0, 4m cos n τ α πL τ n τ K 1 (mL τ n τ ) .
By taking the twist phase as α = 2πa/N and summing over a = 1, · · · , N , we obtain the pseudoentropy for a single scalar field