Properties of the twisted Polyakov loop coupling and the infrared fixed point in the SU(3) gauge theories

We report the nonperturbative behavior of the twisted Polyakov loop (TPL) coupling constant for the SU(3) gauge theories defined by the ratio of Polyakov loop correlators in finite volume with twisted boundary condition. We reveal the vacuum structures and the phase structure for the lattice gauge theory with the twisted boundary condition. Carrying out the numerical simulations, we determine the nonperturbative running coupling constant in this renormalization scheme for the quenched QCD and Nf=12 SU(3) gauge theories. At first, we study the quenched QCD theory using the plaquette gauge action. The TPL coupling constant has a fake fixed point in the confinement phase. We discuss this fake fixed point of the TPL scheme and obtain the nonperturbative running coupling constant in the deconfinement phase, where the magnitude of the Polyakov loop shows the nonzero values. We also investigate the system coupled to fundamental fermions. Since we use the naive staggered fermion with the twisted boundary condition in our simulation, only multiples of 12 are allowed for the number of flavors. According to the perturbative two loop analysis, the Nf=12 SU(3) gauge theory might have a conformal fixed point in the infrared region. However, the recent lattice studies show controversial results for the existence of the fixed point. We point out possible problems in previous works, and present our careful study. Finally, we find the infrared fixed point (IRFP) and discuss the robustness of the nontrivial IRFP of many flavor system under the change of the analysis method. A part of preliminary results was reported in the proceedings [1, 2] and the letter paper [3]. In this paper we include a review of these results and give a final conclusion for the existence of IRFP of SU(3) Nf = 12 massless theory using the updated data.

We report the nonperturbative behavior of the twisted Polyakov loop (TPL) coupling constant for the SU(3) gauge theories defined by the ratio of Polyakov loop correlators in finite volume with twisted boundary condition. We reveal the vacuum structures and the phase structure for the lattice gauge theory with the twisted boundary condition. Carrying out the numerical simulations, we determine the nonperturbative running coupling constant in this renormalization scheme for the quenched QCD and N f = 12 SU(3) gauge theories.
At first, we study the quenched QCD theory using the plaquette gauge action. The TPL coupling constant has a fake fixed point in the confinement phase. We discuss this fake fixed point of the TPL scheme and obtain the nonperturbative running coupling constant in the deconfinement phase, where the magnitude of the Polyakov loop shows the nonzero values.
We also investigate the system coupled to fundamental fermions. Since we use the naive staggered fermion with the twisted boundary condition in our simulation, only multiples of 12 are allowed for the number of flavors. According to the perturbative two loop analysis, the N f = 12 SU(3) gauge theory might have a conformal fixed point in the infrared region. However, the recent lattice studies show controversial results for the existence of the fixed point. We point out possible problems in previous works, and present our careful study. Finally, we find the infrared fixed point (IRFP) and discuss the robustness of the nontrivial IRFP of many flavor system under the change of the analysis method.
A part of preliminary results was reported in the proceedings [1,2] and the letter paper [3]. In this paper we include a review of these results and give a final conclusion for the existence of IRFP of SU(3) N f = 12 massless theory using the updated data.

Introduction
Lattice gauge theory is one of the most successful regularization tools to understand QCD. Since it can be applied even to the strong coupling region, we can investigate the nonperturbative properties of QCD. Based on the success of the lattice QCD, there are several application to the lattice gauge theories with different gauge groups, numbers of flavors and fermion representations. Among recent lattice studies, the search for the conformal or nearly conformal field theory in the infrared (IR) regime has been performed, motivated by both theoretical and phenomenological interests [4] - [48]. If there is an IRFP with nonzero gauge coupling constant, then the low energy physics would show a behavior different from QCD.
The theory is scale invariant but interacting, and the chiral symmetry is preserved in IR region.
In recent studies to search for the infrared fixed point (IRFP), in particular, there are many independent studies in the case of SU(3) gauge theory coupled with N f = 12 fundamental fermions. The existence of the IRFP in N f = 12 theory was predicted by the perturbative beta function at 2-loop [49] and higher [50] in the MS scheme. The phase structure of N f expansion was also studied in the paper [51]. Based on these analytical studies, there are several works using lattice simulations to investigate the nonperturbative running coupling constant and the phase structures, although the results have been controversial hitherto. In Ref. [8], the running coupling constant was computed in the Schrödinger functional (SF) scheme [52][53][54], and exhibited scale independent behavior in the IR at coupling g 2 * SF ∼ 5. And studies with the MCRG method [22,23], studies on the phase structure in the finite temperature system [17,18] and the scaling behavior with mass deformed theory [10,27] show the evidence of the IRFP. In the studies of SU(3) gauge theory with N f = 7 and 10 [4,41], they found that theory is in the conformal window which also suggest that N f = 12 theory is conformal. On the other hand, the studies of the mass scaling behavior [15] and the spectrum of the Dirac operator and the chiral symmetry [12,21] show the evidence that this theory is not conformal at low energy. This situation is confusing, since the existence of the fixed point must be scheme independent. The scheme independence is understood as follows. Let us consider a relationship between two renormalized coupling constant defined in different renormalization schemes, g 1 = F (g 2 ). The beta functions of these two schemes are related as β(g 1 ) = ∂F ∂g 2 β(g 2 ). (1.1) A zero point of the beta function is thus scheme independent except for singular point of the transformations. One possible reason for the controversial situation could be the underestimate of the discretization errors. The nonperturbative running coupling constant using lattice simulations can be obtained using the step scaling method. This method is established by the paper [52] and the nonperturbative running of the renormalized coupling constant in the continuum limit can be obtained. One important point which one should bear in mind is that the careful continuum extrapolation and estimation of the systematic uncertainty are important in the low β (β ≡ 6/g 2 0 where g 0 is the bare coupling constant) region. However, there is no study of the running coupling constant which takes care of the discretization error carefully at least in the case of N f = 12. For example, in the paper [8], the constant continuum extrapolation is taken. That means the discretization effects, which is the renormalization scheme dependent, is neglected.
Another reason may be the bad choice of the value of β . In several previous works, the specific value of β is chosen without any reasons. In the lattice gauge theory with many flavor improved staggered fermion, it was reported that there is a new bulk phase in the strong coupling regime [19,21,24]. Furthermore, the existence of chiral broken phase in the strong coupling limit for the SU (3) theory with N f ≤ 52 is also reported [20]. If the simulation is performed within the bulk phase, it gives an unphysical results because these bulk and chiral 2/41 broken phases are not connected with the continuum limit with asymptotically free ultraviolet fixed point. To avoid such phases, the global parameter search and the determination of the parameter region which is obviously connected to the high β region is needed.
This work reports a study of the phase structure and the running coupling constant for SU (3) gauge theories with N f = 0 and 12. We use the plaquette gauge and the naive staggered fermion actions. Firstly, we study the phase structure of these theories with both analytical and numerical methods, and then compute numerically the running coupling constant with the twisted Polyakov loop (TPL) scheme in the deconfinement phase. The TPL coupling was proposed by de Divitiis et al. [55] for the SU(2) case, and we extend it to the SU(3) theory. This renormalization scheme has no O(a) discretization error, which is of great advantage when we take the continuum limit. Another advantage of this scheme is the absence of zero mode contributions thanks to the twisted boundary condition [56,57]. This regulates the fermion determinant in the massless limit, which enables simulation with massless fermions. In this work, we take the continuum limit carefully, and show the existence of the IRFP in the N f = 12 theory if we include the systematic uncertainty coming from the continuum extrapolation.
This paper is organized as follows. We give a definition of SU(3) TPL renormalized coupling, and the tree level calculation of the TPL coupling in Sec. 2. We show the running coupling constant in the case of quenched QCD theory and the scaling behavior of the scheme and the nonperturbative property of the running coupling constant in Sec. 3. We confirm that the renormalized coupling constant in the TPL scheme approaches to a constant if the theory is in the confinement phase. We discuss how to distinguish such a fake "fixed point" and the true one in this renormalization scheme. In Sec. 4, we discuss the vacuum structure and the center symmetry of SU(N c ) gauge theory, in the system coupled with fermions by the semi-classical analysis in the case of N f = 12 theory to define the true vacua in this setup. We also show the numerical results of the phase structure of massive and massless N f = 12 SU(3) theory to determine the parameter region suitable for IRFP search in Sec. 6. In Sec. 7, we study the running coupling constant for massless N f = 12 case. We firstly study the global behavior of the step scaling function using the data set given in Appendix A, and then investigate the existence of the IRFP by the local fit using the additional data set given in Appendix B in the strong coupling regime. The detailed discussion for the stability of the IRFP is given in Appendix C. We also discuss the taste breaking effects in this analysis in Appendix F.
One of the main results of this paper is that we found the IRFP at and the critical exponent of the β function at the IRFP is There is an independent paper using the similar idea [58]. In the present work, we add the discussion of the quenched QCD and whose fake fixed point in the TPL scheme. Also we give a detailed discussion on the vacuum structure and phase structure in N f = 12 with the twisted boundary conditions and the taste breaking of the staggered fermion, and add the new data showing the strong evidence of the IRFP beyond the systematic uncertainty. Data analysis is also refined in various ways as discussed below. The differences from the 3/41 paper [58], including the discrepancy of the value of g * 2 TPL by more than 2-σ, are discussed in detal in Appendix D and E.

Twisted Polyakov loop (TPL) scheme
One of the nonperturbative definitions for the renormalized coupling constant can be given by a divergence-free ratio (A N P ) of nonperturbative amplitudes. If the tree level value of the quantity is proportional to the squared bare coupling constant as A tree = kg 2 0 , where k is the constant which is calculated by the tree level quantity, then we can define the nonperturbative renormalized coupling constant from the nonperturbative ratio A N P by identifying the renormalization factor of the amplitude as the quantum correction of the coupling constant: Since the lattice simulation gives us the value of A N P , what we have to do is to find a ratio of tree level amplitudes A tree which is proportional to the squared bare coupling constant. Twisted Polyakov loop (TPL) scheme is one of such nonperturbative renormalized coupling schemes defined in finite volume. This scheme is given in Ref. [55] in the case of SU(2) gauge theory, choosing the ratio of Polyakov loop expectation values for twisted and untwisted directions as the quantity A N P . We extend the definition in Ref. [55] to the SU(3) case. Although this scheme can be defined in the continuum finite volume, in this section we start a brief review of the definition of TPL scheme on the lattice.

The definition of TPL scheme in the SU(3) gauge theory
To define the TPL scheme, we introduce twisted boundary condition for the link variables (U µ ) in x and y directions and the ordinary periodic boundary condition in z and t directions on the lattice: for µ = x, y, z, t and ν = x, y. Here, Ω ν (ν = x, y) are the twist matrices which have the following properties: and for a given µ and ν( = µ). The gauge transformation U µ (r) → Λ(r)U µ (r)Λ † (r +μ) and Eq. (2.2) imply Λ(r +νL/a) = Ω ν Λ(r)Ω † ν . (2.4) In the system coupled with fermions, we also have to define the twisted boundary conditions for fermions. Naively, one might think that the twisted boundary condition for lattice 4/41 fundamental fermions can be defined by ψ(x +νL/a) = Ω ν ψ(x), (2.5) for ν = x, y. However, this results in an inconsistency when changing the order of translations, namely, for ρ, ν = x, y. To avoid this difficulty, we introduce a "smell" degrees of fermion N s [60], which can be realized by an integral multiple of the number of color symmetry N c . We identify the fermion field as a N c × N s matrix (ψ a α (x)), where a (a = 1, · · · , N c ) and α (α = 1, · · · , N s ) denote the indices of the color and smell. We can then impose the twisted boundary condition for fermion fields as for ν = x, y directions. Here, the smell index can be considered as a part of "flavor" index, so that the number of flavors should be a multiple of N s , in our case N s should be the multiple of N c = 3. In our simulation, we use staggered fermion which contains four tastes. Now, we can label the flavor(= i, α), where i and α denote the taste of staggered and smell indices respectively. This restricts the number of flavors to multiples of 12 in SU(3) gauge theory with twisted boundary condition. The renormalized coupling in the TPL scheme is defined by taking a ratio of Polyakov loop correlators in the twisted (P x ) and untwisted (P z ) directions: Because of the twisted boundary condition, the definition of Polyakov loops in the twisted directions are modified as, in order to satisfy gauge invariance and translational invariance. At tree level, this ratio of Polyakov loops is proportional to the bare coupling. The factor on the lattice (k latt ) is obtained by analytically calculating the one-gluon-exchange diagram. In our simulation, we choose the explicit form of the twist matrices [61], The Feynman rule for the SU(N c ) gauge theory on the lattice with the twisted boundary condition is given in Appendix B in the paper [55]. The value of k latt is given as whereL = L/a,r = (x, y, z, L/2a) andk µ denotes the momentum in each direction. In the twisted direction,k µ is given by the sum of the physical and the unphysical twisted momenta: where n ph µ = 0, · · ·L/2 − 1 and m ⊥ x,y = 0, 1, · · · , N c − 1 with (m ⊥ x , m ⊥ y ) = (0, 0). The momentumk ⊥ can be identified as the color degree of freedom (N 2 c − 1) in the color basis (see the Appendix B in the paper [55]).
In the continuum limit the proportionality factor (k) for SU(3) gauge theory is calculated analytically as: The values of k latt in Table 1 can be fitted by a linear function of O(a 2 ) instead of O(a), as expected.

The TPL coupling for the quenched QCD
In this section, we study the TPL coupling constant in the quenched QCD. In the first subsection, we discuss the property of TPL coupling in the quenched theory and investigate the parameter region where the study of the TPL coupling makes sense. We obtain the running coupling constant for the quenched QCD in Sec. 3.2.

Phase structure and TPL coupling constant
The TPL coupling constant is defined by taking the ratio of the correlators of Polyakov loop in the twisted and the untwisted directions. If the theory is in the confinement phase the 6/41 correlation length of the Polyakov loop is shorter than the volume, and the gluon does not feel the boundary effect. In such a situation, we can expect that the ratio of the Polyakov loop correlators becomes unity, and give a fake fixed point. For this reason, it is awkward to extract the running coupling and try to give a physical meaning to it in such region. The quenched QCD theory shows the confinement/deconfinement phase transition in the finite volumes, and we can use the TPL running coupling only in the deconfinement phase, where the magnitude of the Polyakov loop shows nonzero values. To see the property of the TPL coupling in both confined and deconfined phases, we study β dependence of the coupling constant at fixed lattice sizes. Apart from discretization errors, the coupling increases as β decreases at a fixed lattice size. In this test, we use smaller lattice sizes, L/a = 2 -6, with relatively low β values. The configurations are generated by the hybrid Monte Carlo algorithm with the Wilson plaquette gauge action. We measure the Polyakov loop and its correlator for every Monte Carlo trajectory, and each data has the same statistics of 20, 000 trajectories.
The TPL coupling constant and the absolute value of the Polyakov loop in t-direction are presented in Fig. 1 1 . The top panels denote the absolute values of the Polyakov loop and the bottom ones denote the corresponding TPL coupling scaled by the coefficient k latt for each lattice size. We found that the absolute value of the Polyakov loop approaches zero in the low energy region. The confinement/deconfinement phase transition occurs at the transition point of β which depends on the lattice sizes. From the bottom panels, we can see the ratio of Polyakov loop (k latt g 2 TPL ) becomes unity below the transition point. Since there can be a fake fixed point due to confinement, there is a question whether we can use this TPL scheme for the conformal fixed point search in IR region. One way to judge that the fixed point is not the fake one is to check the the value of renormalized coupling. Assume that a theory has IRFP. The fake fixed point appears at g 2 TPL ∼ 1/k ∼ 32. If there is an IRFP at g 2 * TPL = 1/k ∼ 32, then we can tell that the fixed point as a physical fixed point. The other important check is to see the phase structure of the theory at the same time. At the true conformal fixed point, the theory must be in the deconfinement phase. There is a possibility of the existence of the bulk phase in the low β region, in which the Polyakov loop shows the confinement and/or chiral symmetry breaking [19,20,24], if the lattice theory contains the dynamical fermions. We discuss this point in the case of N f = 12 in Sec. 6.

Running coupling constant for quenched QCD
Now, we would like to present the results for the running coupling constant in quenched QCD. The main result was already reported in [1]. The gauge configurations in quenched QCD are generated by the pseudo-heatbath algorithm and overrelaxation algorithm mixed in the ratio 1:5. One combination of the pseudo-heatbath and 5 overrelaxation steps called a "sweep" in the following. In order to generate the configurations with the twisted boundary condition we use the trick proposed by Lüscher and Weisz [57]. To reduce large statistical fluctuation of the TPL coupling constant, as reported in Ref. [62], we measure Polyakov loops at every Monte Carlo sweep and perform a jackknife analysis with the bin size of 1 We drop the data of the ratio of Polyakov loop for β = 1.0, L/a = 2, β = 5.0, L/a = 4 and β ≤ 5.5, L/a = 6, since these error bars become huge. They are consistent with 1 as expected.

8/41
To investigate the evolution of the renormalized running coupling, we use the step scaling method [52]. Firstly we choose a value of the renormalized coupling u=g 2 TPL (β, a/L) at the energy scale µ = 1/L. For each L/a in the set of reference lattice size, we find the value of β which produces a given value of the renormalized coupling, u. Then, we measure the step scaling function on the lattice Σ(u, a/L; s)=g 2 TPL (β, a/sL)| g 2 TPL (β,a/L)=u , (3.1) at the tuned value of β for each lattice size sL/a. Here, s is the step-scaling parameter. The step-scaling function in the continuum limit σ(s, u) is obtained by taking the continuum extrapolation of Σ(u, a/L; s): This step scaling function (σ(s, u)) corresponds to the renormalized coupling at the scale µ = 1/sL. To obtain the running coupling constant in the broad range of scale, we identify the value of σ(s, u) with the new input value u at the energy scale µ = 1/sL and repeat the procedures, and then obtain the step scaling function σ(s, u) which corresponds to the renormalized coupling at the lower energy scale µ = 1/s 2 L. Repeating this procedure, we can recursively obtain the renormalized couplings at the scales µ = 1/s n L (n = 0, 1, 2, · · · ). Figure 2 shows the β dependence of the coupling constant in TPL scheme at various lattice sizes. The results are fitted at each fixed lattice size to the interpolating function which is similar to the one used in Ref. [8], where A i are the fit parameters, and 4 ≤ B ≤ 5, n = 3, 4 are employed. As reference lattice sizes of the step scaling, we use L/a = 4, 6, 8, 10. The step scaling parameter is s = 1.5, and we estimate the coupling constant for L/a = 9, 15 from interpolations at the fixed β using the above fit results of all the lattice sizes. We take the continuum limit using a linear function in (a/L) 2 , because the TPL scheme involves no O(a/L) error. We found that the coupling constant in the TPL scheme exhibits a nice scaling behavior even at the smaller lattice sizes, as shown in Fig. 3.
The running of the TPL coupling constant in quenched QCD with 24 steps is shown in Fig.4 together with one-and two-loop perturbative results. The horizontal axis corresponds to the energy scale. The energy scales are normalized at L = L 0 with g 2 (L 0 ) = 0.65. The nonperturbative running coupling constant is consistent with one-and two-loop perturbative results in the high energy region (L 0 /L ≥ 0.1). Comparison with the other nonperturbative analyses in SF scheme and Wilson loop scheme are also interesting. Both the renormalized couplings in SF scheme ( Fig. 1 in the paper [63]) and Wilson loop scheme ( Fig. 8 in the paper [64]) show similar behavior; i.e., it runs faster than the one in one-loop result in the nonperturbative region. On the other hand the TPL running runs slightly slower than the one-loop result.
From this quenched study, we conclude that we can control both the the statistical and systematic errors of the TPL coupling constant. Furthermore we find that the TPL coupling constant in quenched QCD has a robust scaling behavior even in a small lattice size, which was also observed in the previous quenched SU(2) calculations [55,62] 4. Vacuum structure of the N f = 12 theory with twisted boundary condition Next, we consider the SU(3) theory coupled to N f = 12 fundamental fermions. In this section, we discuss the center symmetry of SU(3) gauge group to define the true vacuum in this theory. The generators of the center symmetry of SU(N c ) pure gauge theory are z = exp(2πil/N c ), where l = 0, 1, · · · , N c − 1. This symmetry is broken by adding fermions to the theory, leading to the existence of true vacua in an SU(3) gauge theory involving massless fermions in the deconfinement phase. We discuss the vacuum structure which is important to the study of the gauge theories in finite volume.
Let us focus on the center symmetry of this theory. Although the Wilson gauge action is invariant under the following transformation for the link variable for each direction, the fermion is not invariant. At the perturbative one-loop level, the semi-classical free energy for the gauge configuration {U } is given as With the twisted boundary condition, the flat potential due to the toron contribution is lifted because of the nonzero momenta in twisted directions and the free energy has 3 4 = 81−fold degenerate classical minima at U µ = exp(2πiθ µ /3)I, where θ µ = 0, 1, 2 for each direction.
The Wilson gauge action (S g ) and the one-loop contribution from gauge part (S one-loop g [U ]) respects the Z 3 symmetry, so that we do not consider them in what follows.
Let us consider the fermion determinant. In the momentum space, there are the physical and unphysical momenta (k µ =k ph µ +k ⊥ µ ) in the twisted directions, that also appear in the gauge field momenta (Eq. (2.12)). In the case of the fermion field, the color and smell degree of freedom of ψ a α in Eq. (2.7) can be transferred into the unphysical momentum degrees of freedom: their number is N c × N s 2 . Here, we replace the momentum aŝ The Z 3 transformation in Eq. (4.1) can be defined on each lattice site independently, so that we can take a typical gauge in which U µ = exp(2πiθ µ /3L)I for whole lattice volume. Then the fermion action in the vacuum U µ = exp(2πiθ µ /3L)I is obtained by the above replacement. The fermion determinant in finite volumeL 4 is thus In Table 2, we give the results of the fermion determinant (Eq. (4.4)) for L/a = 6. We find that the vacuum free energy is independent of θ x,y and there are three types of vacua classified with θ z,t . The first one is a "trivial vacuum", in which vacuum (θ z , θ t ) = (0, 0). This vacuum has 9−fold degeneracies. The value of the free energy is highest among three types of vacuum, so that it will decay to the true vacuum. The second one is a "half-trivial vacuum", in which one of θ z,t is 1 or 2 and the other one is θ µ = 0. This vacuum has 36−fold degeneracies. The free energy is higher than the one of the third vacuum, so that this vacuum is also unstable. The third one, in which the free energy is lowest, is a "non-trivial vacuum". Both θ z and θ t take 1 or 2, and there are also 36−fold degeneracies. The free energy has minima at this vacuum where all classical link variables for z and t directions has a nontrivial phase U z,t ∝ exp(±2πi/3L). This means that the Polyakov loop for z direction also has a non-trivial phase exp(±2πi/3).
This classification holds for generic lattice size, and it turns out that the difference of the free energy between the true non-trivial vacua and the other vacua becomes small as the lattice size becomes larger, where as the potential barrier becomes higher. If we change the fermion boundary condition in z and t directions, the semi-classical free energy has minima at other vacua.

Simulation setup for the N f = 12 theory
Our numerical simulation is performed in the following setup. The gauge configurations are generated by the hybrid Monte Carlo algorithm with the Wilson gauge and the naive staggered fermion actions. In Sec. 6, the simulations are carried out with lattice sizes L/a = 4, 8 and 12 at several low β and a broad range of ma to study the phase structure in this system. We generate 1, 000 -2, 000 trajectories for each parameter set in the case of L/a = 4 and 8, and also generate 500 -1, 000 trajectories for each in the case of L/a = 12.
The measurement of the coupling constant in Sec. 7 are carried out with lattice sizes L/a = 6, 8, 10, 12, 16 and 20 3 at around thirty values of β in the range 4.0 ≤ β ≤ 100 with ma = 0. To reduce statistical fluctuations, we measure the Polyakov loops at every trajectory and bin the data by taking the autocorrelation into account. Using the jackknife method, typical statistical errors of correlator are found to be 2 -3%. We also estimate the statistical error within the bootstrap method in the whole analysis in Sec. 7, and obtained consistent results.
6. Simulation results: Phase structure of N f = 12 SU(3) theory with the twisted boundary condition In this section, we investigate the phase structure of the N f = 12 fermion theory on the lattice. Although the main purpose of this paper is to study the running coupling in the massless theory, in order to fully understand the phase structure we need to understand the phase structure of the whole region of the theory space including the mass parameter. Actually, there are several studies which reported the existence and absence of the bulk phases in the case of N f = 12 staggered fermion system [19,24]. In the paper [24], the authors found the existence of a bulk phase where shift symmetry and chiral symmetry are weakly broken, and the paper [19] reported that it is caused by the improvement of the staggered fermion. Furthermore, there is the spontaneous chiral symmetry broken phase in the strong coupling limit for N f ≤ 52 [20,65] for SU(3) massless fermion theory. In our simulation, we use the naive staggered fermion and introduce the twisted boundary condition.
It is important to show the phase structure in our lattice setup independently to identify the parameter range suitable for the study the TPL coupling constant. In Sec. 6.1 and Sec. 6.2, we study the β and volume dependence of the expectation value of the plaquette and the Polyakov loops from which we determine the phase structure of the theory in the β − ma plane. In Sec. 6.3, we discuss a contribution of the vacuum tunneling due to the lattice artifact to the TPL coupling constant. Due to the finite lattice spacing, the tunneling behavior between these vacua occurs during Monte Carlo simulation and it must be a lattice artifact. As we have shown using semi-classical analysis in Sec. 4, there are the degenerate vacua in this lattice setup. We show that the contribution is negligible in our simulation.   First of all, we found there are discontinuities in the plaquette as a function of β along the lines of fixed ma. The discontinuity appears between β = 4.8 and β = 4.6 for ma = 0.2 and in the smaller mass region it appears at the lower β. Near the massless limit the gap exists around 3.6 ≤ β ≤ 4.0. 13/41 The small (red) arrows near the massless at β = 3.8 on the left panel in Fig. 5 shows the detailed histories of the thermalization. We find that there are two different values in the 0 ≤ am ≤ 0.0125 region. The configurations giving larger values of the plaquette at the same ma are generated starting from configuration with massless fermions in β = 4.0; on the other hand those giving smaller values are obtained starting from the configuration with massive fermions at fixed β. The hysterisis clearly indicates that there is a first order phase transition around this region. At β = 4.0 and β = 3.6, there is no dependence on the thermalization process.
We also study larger mass region. The right panel in Fig. 5 is the same plot as the left one for a broader region of ma. In the quenched limit, we know that there is the first order phase transition. In the figure, we plot the data for the quenched lattice at ma = 1.0. The gap seems milder in the larger mass region.  We also investigate this first order phase transition near the massless region by changing the lattice volume in Fig. 6. There are slight differences between L/a = 4 and L/a = 8 for the critical value of β. For example at ma = 0.175 the data for β = 4.6 on (L/a) 4 = 4 4 is in the upper phase of the gap (See Fig. 5), but it moves to the lower phase in the case of (L/a) 4 = 8 4 (See left panel in Fig. 6). The results on (L/a) 4 = 8 4 and 12 4 show the similar behavior at least the present interval of β and ma (∆β = 0.2, ∆ma = 0.01 -0.025), and there is no clear volume dependence. Since the massless simulation at β = 3.8 needs extremely finer 14/41 molecular-dynamics time step size than ∆τ = 0.002 (τ = 1 is 1 trajectory), practically we could not generate the data. The position of β where the simulation becomes quite costly is the same for both L/a = 8 and L/a = 12. It suggests that near the massless region there is a bulk phase transition in β < 4.0.

Polyakov loop
Next, let us investigate the Polyakov loop. Since the dynamical fermions breaks the center symmetry explicitly, there is no clear order parameter for the deconfinement phase transition. However, here we use the word "deconfinement" phase for the region in the theory space where magnitude of Polyakov loop is clearly nonzero on the lattice.
In the case of the quenched QCD, there is Z 3 symmetry, and there is tunneling between them on finite lattices. In the case of N f = 12 massless theory with the twisted boundary condition, the Z 3 symmetry is broken and the true vacua is the one that the Polyakov loops in the untwisted directions have the nontrivial phase as explained in Sec. 4.
At first, we present the scatter plots of the typical Polyakov loops for the massless theory in the twisted and untwisted directions in the left and right panels in Fig. 7 respectively. The we confirm that the Polyakov loops in both directions are consistent with the results of from the semi-classical analysis in Sec. 4 even in the strong coupling region. The effect of the tunneling on the TPL coupling will be discussed in the next subsection. Clearly there is a gap for each ma. In the case of the massless fermions, the measured data at β = 3.8 are shifted; one of them shows at β = 3.77 which corresponds to the configurations at β = 3.8, ma = 0 whose plaquette value is the larger one in Fig. 5. The other shows at β = 3.83 which corresponds to the configurations at the same β and ma whose plaquette value is the smaller one in Fig. 5. Right panel: the same plot in the case of (L/a) 4 = 8 4 .
Next, we show the real part of Polaykov loop in t-direction. The left two panels in Fig. 8 show the ones at several fixed ma and β in the case of L/a = 4. We can find that there is a gap of the real part of the Polyakov loop at fixed ma data, and the value of β at the gap corresponds to the critical value of β of confinement/ deconfinement phase transition. In the case of massless fermions on L/a = 4 we find a gap at the β = 3.8. For β smaller than the gap position the real part of the Polyakov loop is not consistent with zero, but it goes to zero continuously. In the finite mass region, there is a weak jump, and the gap become larger in the smaller mass region. The value of the critical β in which the data shows the jump is the same with the plaquette study.
Let us study the lattice size dependence. We show the real part of the Polyakov loop in the case of L/a = 8 in the right panels in Fig. 8. There is no clear jump in the case of L/a = 8, but the real part of the Polyakov loop approach to the zero in the low β region. We define the critical value of β as the largest value of β for which the real part of Polyakov loop becomes consistent with zero. Again, the critical value of β is the same with the plaquette study. We can conclude that in Figs. 5 and 6, the phase for β larger than the gap position can be identified as the deconfinement phase and that for β smaller than the gap position is the confinement phase.
Note that there is one misleading exceptional data at β = 6.0, ma = 0.4, L/a = 8. The real part of Polyakov loop is consistent with zero. However, we can find that the theory is in the 16/41 deconfinement phase from the scatter plot of the Polyakov loop on the complex plane. Since at the parameter the fermion mass is too heavy, there is the tunneling between Z 3 vacua as in the quenched case. That is the lattice artifact, so that the data does not say the theory is in the confinement phase.
From the comparison of the data at L/a = 4 and L/a = 8, we found there is a volume dependence of the critical value of β in the massive region, while there is no dependence near massless region. Actually, in the top panels of Figs. 8, the theory with ma = 0.4, L/a = 8 is in the confinement phase for β ≤ 5.0, while in the case of L/a = 4 the theory is in the deconfinement phase for 4.8 ≤ β. Figure 9 shows the scatter plots of the Polyakov loop on the complex plane. These plots show that there is a clear volume difference. On the other hand, in the case of the (nearly) massless fermion, the transition point does not show the volume dependence, indicating that the transition is bulk one. The theory with the massless fermion in both L/a = 4 and 8 lies in the deconfinement phase for β ≥ 4.0. At least the current interval of β and ma, we cannot find the volume dependence in the small mass region. Since in the quenched limit we know that there is a finite volume phase transition, we expect that there is a crossover for both bulk and the finite volume phase transition in the middle range of the fermion masses.
Furthermore, in the case of L/a = 12, we cannot find clear nonzero value of the Polyakov loop in the whole region since the current preliminary statistics is small and the Polyakov loop in the low β region is noisy. The gap of the plaquette and the confinement/deconfinement phase transition seems to occur simultaneously, and there is no difference between L/a = 8 and L/a = 12 on the plaquette behavior. At least we found that the position of β where the simulation becomes quite heavy is the same in the case of L/a = 8.
Finally, we study the phase structure of the massless fermion N f = 12 QCD for β ≥ 4.0 and L/a = 6-20. We find that all configurations, which are used for the running coupling constant study in Sec. 7, live in the deconfinement phase (although it might be trivial since the transition seems to be the bulk and we concentrate on the parameter region within β ≥ 4.0). Figure 10 shows the real part of the Polyakov loop for the low β region for each 17/41  The summary of the phase structure and the available region of the TPL coupling for the quenched and the massless N f = 12 QCD is the following. Figure 11 is a sketch of the phase structure for the naive staggered N f = 12 SU(3) theory. In the case of the quenched QCD, the correlation length becomes shorter in the lower β region, and there is the finite volume phase transition where the theory goes to the confinement phase. In the case of the massless N f = 12 SU(3) theory, there is the similar behavior while the transition seems to be the bulk 18/41 one at β < 4.0. In the study on the running coupling constant in TPL scheme, we should focus on only the deconfinement phase on the lattice. Furthermore, in β ≥ 4.0 region with massless fermions, we also investigate the eigenvalue of Dirac operator which is presented in Appendix F. The lowest eigenvalues are clearly nonzero even in the lowest β for all lattice sizes. According to the Banks-Casher relation, it implies that the chiral symmetry is restored in β ≥ 4.0. In Sec. 7, we finally find the IRFP at higher β values than the bulk phase transition point, although the values of physical critical β at physical IRFP depend on the lattice sizes. Our phase diagram (Fig. 11) is completely consistent with the conjectured phase diagram in the paper [20] (Fig. 10) by Ph. de Forcrand et al., where they study the strong coupling limit.

Tunneling behavior between true vacua
During the Monte Carlo simulation, the tunneling can occur between the two degenerate vacua in each untwisted direction independently. The tunneling behavior is a lattice artifact, since the potential barrier between the vacua becomes finite at the finite lattice spacing. In this subsection, we consider how the TPL coupling is disturbed by the tunneling behavior.
The tunneling is expected to occur more frequently in the strong coupling region. We observe some tunnelings in low β region, although the number of the tunneling is quite small, and decreases as β increases. For L/a = 6 we observe the tunneling 7 times per 6, 000 trajectories at β = 4.0, and only 3 times per 90, 000 trajectories at β = 6.0.
A typical example of the tunneling is shown in figure 12. The left panel denotes the history of the imaginary part of the Polyakov loop in the untwisted direction at β = 4.0 and L/a = 6. As can be seen in Fig. 12, the sign of the imaginary part is flipped at around the 1, 400 -1, 600th trajectories.
Let us see the effects of the tunneling on the TPL coupling defined by the ratio of the Polyakov loop correlators in Eq. (2.8). The right panels in Fig. 12 present the histories of the numerator and denominator of the TPL coupling obtained from the same configuration as in the left panel. During the tunneling, i.e., the 1, 400 -1, 600th trajectories, both the results seem not to exceed the range of each statistical fluctuation. This means that the tunneling does not significantly affect the result of the TPL coupling calculated by the ratio of the expectation values for the numerator and denominator.
Therefore we conclude that effects of the tunneling on the TPL coupling is negligible at least in our calculation due to the rare occurrence of the tunneling and the large statistical fluctuation of the Polyakov loop correlators.

Simulation results: Step scaling function for N f = 12 SU(3) gauge theory
In this section, the step scaling function in the case of N f = 12 flavor is derived. As in the quenched QCD case, we use the step scaling method to find the IRFP. In this section, we focus on another quantity, which is called the growth rate of step scaling function, rather than the running coupling constant. As we explained in Sec. 3.2, the running coupling constant is derived using "recursively" the step scaling procedure. Since the error of σ(s, u) feeds in to the input renormalized coupling (u) in the next step scaling procedure, the errors from σ(s, u) accumulate in the running coupling towards lower energy. On the other hand, the step scaling function σ(s, u) with no error accumulation can be defined independently for 19/41  each step and the growth rate σ(s, u)/u becomes unity when there is a zero in the beta function. Therefore the growth rate σ(s, u)/u is a suitable quantity for searching the IRFP. At first, we will discuss the global behavior of the growth rate from the perturbative to the IR region in Sec. 7.1. The nonperturbative running behavior shows the signal of the conformal fixed point in the IR region. Then, we focus on the low energy region only and derive again the step scaling function by using the data only in the strong coupling region in Sec. 7.2. We discuss the stability of the IR fixed point by considering several systematic uncertainties and derive the universal quantity for the exponent of the beta function around the IRFP in Sec. 7.3.

The global behavior of the step scaling function
Before explaining the simulation details, we address the guiding principles of our simulation to show the global behavior of the running coupling. We use the step scaling method as in the quenched case. For each L/a we interpolate the data in β. Since the choice of the interpolation function can affect the step scaling function, we should take care of the following two points: (1) We generate data in a broad range of β, with intervals such that the renormalized coupling constant (g 2 R ) grows almost constantly in each interval. Thus the interval of β is large in high β region while small in the low β region. Each data has a similar accuracy (2 − 3%).
(2) We employ fit functions for β interpolation which reproduce the tree level result g 2 R ≃ g 2 0 on each lattice size in extremely high β region. These guiding principles ensures the stability of our fit results under the change of fit functions and the number of data. Point 1 is needed to ensure that the fit functions do not favor any special region of the data when we interpolate our data in β or extrapolate to in (a/L) 2 . To satisfy this point is very important to search for the IRFP by using the numerical simulations, since finally the interpolation function of the data decides the step scaling function. We discuss the dependence of the data set for the global fit analysis in Appendix D, in the case where the data are concentrated in some particular region. The result in Appendix D do not have a nice agreement with the perturbative result and in the IR region the position of the IRFP strongly depends on the fit range. Point 2 is needed to reduce the effect of statistical fluctuation in high β region. In high energy region, since the coupling runs very slowly, we need extreme high statistics to reproduce the perturbative result only by the numerical data. The assumption of the point 2 makes the result stable in high energy region. In Fig.13, we show our simulation results for the renormalized coupling in TPL scheme as a function of 1/β for each L/a. The raw data are given in Tabs. A1 and A2 in Appendix A. The left panel in the Fig. 13 shows a global behavior of the TPL coupling. We can see the high energy behavior seems to be almost linear in 1/β as expected from the perturbation theory. The right panel focuses on the low β region. In the low β region for L/a = 6 the TPL coupling has a maximum at β = 4.3. In contrast to the Schrödinger functional scheme [8], the renormalized coupling gets larger for larger volume for all the range of β. We consider that this difference comes from the lattice artifact which depends on the renormalization scheme. To remove the effect, the careful continuum extrapolation is necessary.
For β-interpolation, we use the following form of fitting function: where C i (a/L) are the fitting parameters and N is the degree of the polynomial. Here, N = 3 -5 are employed. We drop the data at β = 4.0, L/a = 6 from the fit to avoid the double solutions when we solve the β values to reproduce the input renormalized coupling (u). This limits us to study the step scaling in the range u ≤ 2.94, where u = 2.94 is the renormalized coupling constant g 2 TPL at β = 4.3, L/a = 6, in this lattice set up. 21/41 To investigate the evolution of the renormalized running coupling, we use the step scaling method as explained in Sec. 3.2. In this study, we take s = 1.5, and denote σ(u) ≡ σ(s=1. 5, u) in the rest of this paper for simplicity. The set of small lattices is taken to be L/a = 6, 8, 10, 12, therefore, we need values of g 2 TPL for L/a = 9, 12, 15, 18 to take the continuum limit in Eq. (3.2). For L/a = 9, 15 and 18, we estimate values of g 2 TPL for a given β by the linear interpolation in (a/L) using the data on the lattices L/a = {8, 10}, {12, 16} and {16, 20}, respectively. To estimate the systematic error of these interpolations, we also performed the linear interpolation in (a/L) 2 , and found that the difference in the results with interpolations in a/L and (a/L) 2 is negligible. Furthermore, we will also show the step scaling with s = 2 in which there is no interpolation of L/a, and the difference of the result should give an indirect estimation of the systematic error from the interpolation in L/a.  In Fig. 14, we show the examples of the continuum extrapolation for obtaining σ(u) in the weak, intermediate and strong coupling regions. We determine the central value of σ(u) by the linear extrapolation in (a/L) 2 with four points; L/a = 6, 8, 10, 12 since the leading discretization error is of O(a 2 ) in this scheme. Note that, in the strong coupling region, each lattice data Σ(u, a/L; s = 1.5) (black data point) is quite larger than u, however, in the continuum limit, σ(u) gets close to u. This indicates that it is very important to take the continuum limit carefully in this kind of analysis. We explain the reason why this continuum extrapolation is chosen as the best in the analyses in Appendix E. We also discuss the taste breaking in the continuum limit in Appendix F. Now, we obtain the step scaling function explained above in a wide range of u. Figure 15 shows the growth rate of the renormalized coupling (σ(u)/u) as a function of u with statistical error which is estimated by jackknife method. We also carried out the bootstrap analysis independently, and found that the results are consistent with this jackknife analysis.
We found two things from this plot. The first one is that the result is consistent with perturbation theory in the weak coupling regime. The TPL coupling coupling with this 22/41  Table 3 The fit range of the local fit analysis for each lattice size. The data L = 9, 15 and 18 are obtained by the L/a interpolation at the fixed β value as explained in Sec. 7.1.
lattice set up looks promising under this analysis method. The other one is the central value of σ(u)/u becomes unity around u = 2.7, demonstrating the signal of a fixed point. This is the first zero of the beta function from the asymptotically free regime. It suggests the existence of an infrared fixed point around the region. Unfortunately, the upper values of the error bars do not cross the line σ(u)/u = 1. This means that we cannot exclude the possibility for the coupling constant to continue growing within the error bar. We will investigate this quantity again by focusing only the strong coupling region and adding the data. Furthermore, will give an estimation of the systematic error of the fixed point coupling in the next subsection.

Low energy behavior and stability of the IR fixed point
In the previous subsection, we found a signal of the IRFP around u = 2.7 from the global fit of the data. Now we focus on the strong coupling region and will determine the fixed point coupling and the related universal quantity. In this subsection, we take a narrow β range in which β-dependence of g 2 TPL can be approximated by linear or quadratic functions of β. We add more data, which is a part of the data shown in Appendix B, to obtain the precise result and discuss the systematic uncertainties of the IRFP.
Practically, we will carry out the step scaling again with the data only in low β region u ≥ 2.0. This region roughly corresponds to the range β ≤ 7.0. We choose the fit range for each lattice size as shown in Table 3. The fitting function is chosen as a simple unconstraint 23/41 polynomial function: where N is the degree of the polynomial and here we take N = 3 for L/a = 6, 8 and N = 2 for the other lattice size. We derived the step scaling function by using the same procedure as in the previous subsection. The growth rate of the step scaling function is shown in Fig. 16.
As a central analysis with solid blue error bar, we take the four point linear extrapolation in (a/L) 2 with statistical error estimated by the jackknife method. The dot (black) error bar includes the systematic error, which we will discuss later. This local fit result clearly crosses the line σ(u)/u = 1, which shows the existence of the IRFP. Figure 17 shows the comparison between the local fit result with the global fit result. These two central value are consistent with each other within 1-σ, despite the change of the data set, the fit range, and the fitting function. The results strongly consolidate the existence of the stable fixed point around u = 2.7. The error bar for the local fit analysis is smaller owning to the additional precise data. We also report the data set dependence and the fit range dependence independently in Appendix C. Now, we would like to estimate the systematic error in our analysis. There are two possible dominant sources of the systematic error. One is from the choice of the fit range for the βinterpolation (Eq. (7.2)). As shown in Fig. 17, there is a small difference between the global fit and the local fit in Table 3, and we also investigate narrower range of the β. Even if we take only the data of β ≤ 7.0 for all lattice size, the fitting function does not show a large difference. We can conclude the systematic error from the choice of the fit range is small in this analysis (see Fig. C1 in Appendix C).
The other dominant systematic error comes from the continuum extrapolation. In Fig. 18, we show the comparisons of several types of continuum extrapolation for u = 2.0, 2.686 and 2.85. As the central value, we take the linear extrapolation in (a/L) 2 for L/a = 6, 8, 10, 12. We 24/41  respectively. The red line shows the extrapolation function linear in (a/L) 2 for 3 data points without the coarsest lattice data. In the case of u = 2.0, the step scaling function is larger than the input value, however, it becomes consistent with u at u = 2.686 and for the larger u it is smaller than the input renormalized coupling constant.
estimate the systematic error by taking the difference between the central value and the result from linear extrapolation without the data on the coarsest lattice L/a = 6. Furthermore we compare the central value with the quadratic extrapolation with all the data at four values of L/a. Figure 18 shows the TPL renormalized coupling has a small systematic error in the 25/41 strong coupling region, and all the values in the continuum limit agree within 1-σ statistical errors. The total error in Fig. 16 is estimated by adding the difference between the continuum extrapolations as a systematic error to the statistical error in quadrature. We conclude that the existence of the IRFP is stable in this analysis. Finally the renormalized coupling at the IRPF is u * = 2.69 ± 0.14 (stat.) +0 −0.16 (syst.).
Here, the jackknife error of the running coupling is used as a statistical error and we estimated the systematic error coming from the continuum extrapolation. The corresponding β value for each lattice size at u * = 2.69 can be calculated from the β interpolation function in Eq. These are the parameter sets which reproduce conformal physics after taking the continuum limit.
In addition, we mention the numerical stability of our results. For this purpose, we perform another step-scaling analysis based on s = 2 with L/a = 6, 8, 10. The continuum limit is taken by linearly extrapolating three points in (a/L) 2 . The growth rate of the renormalized coupling is shown in Fig. C1 in Appendix C. We find that the running behavior in s = 2 step scaling is also consistent with that in the case of s = 1.5, and the IRFP is found at u * = 2.49 ± 0.19 (stat.). This consolidates the existence of the IRFP in this theory.

Critical exponent
Finally we will derive the critical exponent at the IRFP. In this theory, we have one irrelevant parameter, which is the renormalized coupling constant, around the nontrivial fixed point. The value of the renormalized coupling at the fixed point depends on the renormalization scheme. If we denote the scheme transformation from one renormalization scheme u 1 to some other one u 2 = F (u 1 ), then in the perturbative region, the function F (u 1 ) can be expanded as a polynomial function. The beta function for the renormalized coupling is universal up to two-loop order, however, the nonperturbatively the one for u 2 scheme is related with the other one as follows: In the vicinity of the IRFP, the beta function can be approximated by Although the value of renormalized coupling at the IRFP is scheme dependent, we can easily find the coefficient γ * g is the scheme independent quantity using Eq. (7.4). Now, we compute γ * g from the slope of σ(u)/u against u, and obtain s −γ * g = 0.79 ± 0.11(stat.) in the central analysis in the Fig. 16. This leads to γ * g = 0.57 +0.35 −0.31 (stat.) +0 −0.16 (syst.), (7.6) where the first error is statistical error using the jackknife method and the second one is the systematic error from the continuum extrapolation estimated by the comparison to the 3 point linear continuum extrapolation. The value of γ * g is sensitive to the variation of the slope, which causes rather large statistical error. For the s = 2 step scaling, the critical 26/41 exponent of the beta function can be derived γ * g = 0.31 +0.21 −0.18 (stat.). This is also consistent with our main results with s = 1.5.
Our result is consistent with γ * 2−loop g ∼ 0.36 and γ * 4−loop(MS) g ∼ 0.28 as estimated using 2-loop and 4-loop (MS scheme) perturbation theory [50,66,67]. The result in the SF scheme is γ * SF g = 0.13 ± 0.03 [8]. Our result provided a value larger than that in the SF scheme. We can conclude that both results are almost consistent with each other since the discrepancy of γ * g is slightly larger than 1-σ. Another scheme independent quantity which is interesting to observe is the mass anomalous dimension at the IRFP. That is the critical exponent for the relevant operator around the IRFP. We will report it in a forthcoming paper [68].

Summary
We gave the explicit definition of the TPL renormalized coupling for SU(3) gauge theory and studied the running coupling constant in the case of quenched QCD and N f = 12 theories. The definition is the extension of the SU(2) gauge theory, and we provided the perturbative calculation to define the coefficient in the case of SU(3).
Firstly we show the TPL running coupling constant in the case of quenched QCD. In the theory, there is confinement/deconfinement phase transition because of the finite volume effect and we study the behavior of the TPL renormalized coupling constant in the both phases. The TPL scheme has the remarkable property that in the extremely low energy limit the coupling constant approaches to the constant (g 2 TPL ∼ 32 in the case of SU(3) gauge theory), when the theory is in the confinement phase. From this analysis, the TPL coupling is found to be useful only in the deconfinement phase, so that we should study the phase structure in the parameter spaces and search for the available region before the running coupling study. The running coupling constant is consistent with the perturbative result in high energy region, and it runs more slowly than that in the 1-loop perturbation in the low energy region. The running coupling constant in nonperturbative region is scheme dependent and is different from SF and Wilson loop schemes.
In the case of SU(3) gauge theory with 12 massless Dirac fermions in the fundamental representation, we inverstigated the phase structure on β-ma space and the vacuum structure related with Z 3 center symmetry to identify the true vacua. We revealed the phase structure for N f = 12 massive and massless fermion theories and found that there is a bulk phase transition near the massless region at a point of β < 4.0 region. In such phase, the TPL coupling is not available since the theory shows the confinement behavior. We used the configuration only in the deconfinement phase to investigate the running coupling constant, and also found that the chiral symmetry seems to be preserved.
We also discussed the vacuum structure and the center symmetry breaking in our simulation setup using the semi-classical analysis, and generated the configurations in the true vacua. The vacuum structure depends on the boundary condition of the fermions, in the case of our definition, the configurations whose Polyakov loop in the untwisted direction has the nontrivial phase shows the minimum of the potential.
Finally, we have found a solid evidence for the existence of an IRFP using the TPL scheme. The coupling constant at the IRFP is g * 2 TPL ∼ 2.69. The stability of the fixed point 27/41 is discussed, and we can conclude there is the IRFP after the systematic uncertainties are included. C. Data set dependence, fit range dependence and the step scaling size dependence We show the two kinds of result for the growth rate from the global fit analysis and local fit analysis in Sec. 7.1 and 7.2 respectively. They are consistent with each other within 1-σ although they have the difference data sets, the different fit range and fitting function each other. Here we would like to show the each systematic uncertainties. The left top panel in Fig. C1 shows the data set dependence. The central analysis includes both Data sets in Appendix A and B, while the blue result is obtained by only Data sets in Appendix A, which is the same data sets with the global fit analysis. The results are consistent with each other within 1-σ. The right top panel in the Fig. C1 shows the fit range dependence. The blue result is obtained by the data in the narrow β range; β ≤ 7.0 for all lattice size. We can find the result is completely consistent with each other, and this local fit analysis is quite stable under the change of the fit range. Actually, we focus on the local β region where all data can be fitted by the linear or quadratic function in β. Such results can be expected when the data can be fitted well.
The bottom panel in Fig. C1 shows the result of s = 2 step scaling. Although the step scaling function depends on the step scaling parameter, the comparison for the position u * must be independent of the step scaling parameter. We can find the position is consistent with that for s = 1.5, so that we can confidently conclude the existence of the fixed point and the interpolation works well in the s = 1.5 step scaling. 32

D. The effect of the nonuniform data on the global fit
The fixed point in this paper shows u * = 2.69, however, the paper [58] shows u * ∼ 2.0. We consider that this difference comes from the mismatching the analysis method and the data sets quality in the paper [58]. In this appendix, we would like to consider the problem. The data in the paper [58] are strongly concentrated around β = 5.0 and they are a part of Data set A and all data of set B. In the paper [58], the authors carried out the global fit analysis. As we mentioned in our guiding principle point 1, when we use the global fit in the broad β region by using a single fitting function, ideally the data do not favor a specific region. Our analysis for the global fit used only the data A in which each data has a similar accuracy and the interval of β is chosen to realize the almost constant growth rate of the renormalized coupling on the lattice. On the other hand, the data set B is strongly concentrated around β = 5.0 -6.0 and a part of the data has quite high statistics in the high β region and for the small L/a. In this appendix, we would like to derive the global step scaling function by using the data in the both sets A and B and discuss the effect of the prejudiced data.
In Fig. D1, the results by using the data set A and the ones by using both A and B are shown in the top and bottom panels respectively. Each procedure of the step scaling is the same with the Sec. 7.1 and Sec. 7.2 for the global and local fit analysis respectively. Each red result in fig. D1 denotes the global fit analysis and the blue one denotes the local fit 33 analysis. The latter global fit result looks very similar with the result of the paper [58]. It shows the worse matching with the perturbative result, although still it is consistent with the perturbative prediction within the statistical error. The values u * becomes smaller than the other ones with more than 2-σ discrepancy. Now, we should consider which result is reliable. The global fit with the broad regime is always dangerous since the interpolated value includes non-vanishing contributions coming from the far region. To remove such effect we also carried out the local fit only with the focusing regime. The comparison the global fit result and the local fit result with both data A and B shows the larger fit range dependence rather than our central analysis in Fig. 17. The guiding principle point 1 must be important for such global fit analysis. That is the reason why we did not use the whole data of the data set B, when we would like to show the global behavior of the TPL running coupling constant.

E. Comments on the estimation method of the discretization effect
The discretization effect of the step scaling function Σ(u, a/L; s) has two origins. The first one is a simple discretization effect of the renormalized coupling in the larger lattice size (sL/a). The second one comes from the tuned value of β to reproduce the input quantity u in the smaller lattice (L/a). When we take the continuum limit, we fixed the physical box size L and the lattice spacing a (= β), and then the leading term of the former (O((a/sL) 2 )) 34/41 is smaller than the later one (O((a/L) 2 )). So, it must be safe to avoid the interpolation of the small lattice sets if we consider the discretization effect seriously.
In our simulation, we have L/a = 6, 8, 10, 12, 16 and 20. Then in the s = 1.5 step scaling we can take 4 data points to estimate the O(a 2 ) effect with avoiding the small L/a interpolation, on the other hand s = 2 step scaling has only 3 data points. One of the advantages to use s = 1.5 step scaling is that there is the finest lattice data(L/a = 12). Furthermore the chi square fit with only one degree of freedom is strongly disturbed the statistical fluctuation, then we take the 4 points linear extrapolation as the central analysis in this paper.
In the paper [58], they carried out the interpolation for L/a = 7 by using L/a = 6 and 8. However, the interpolation is dangerous since it is the coarsest two lattices interpolation. Actually, the raw data of the TPL (Figs. 13) shows the largest difference between L/a = 6 and 8 in the low β region and that must induce a large uncertainty of the interpolation. Furthermore, the estimation of the systematic uncertainty between 4 data points linear extrapolation for L/a = 6, 7, 8, 10 → 12, 14, 16, 20 and 3 data points linear extrapolation for L/a = 7, 8, 10 → 14, 16, 20 is nonsense, since the data of L/a = 7 is generated by L/a = 6 lattice data and thus the later extrapolation does not remove the effects of the coarsest lattice.

F. Eigenvalue of the Dirac operator
In this appendix, we report results for the eigenvalue distribution of the Dirac operator. We confirm two things from the quantity. At first, we show the global shape of the eigenvalue distribution. We find that the data in the weak coupling region, where the perturbative theory gives good approximations, is consistent with the tree level analysis in Sec. F.1. Furthermore, the β dependence of the data is smooth in the whole region and the lowest eigenvalue is nonzero even in the lowest β in our simulation parameter. These observations also indicate that the theory is in the deconfinement and chiral symmetric phase. Secondly, we discuss the taste breaking in our simulation in Sec. F.2. In the case of the high β and the large lattice extent, the raw data of the low lying eigenvalues shows the degeneracy of the taste. In the strong coupling region, we find that the effect of the taste breaking becomes mild in the continuum limit. We consider up to the order a 4 for the continuum extrapolation, which is the same order as for the running coupling study. How much the taste breaking recovers up to this order would be an indirect estimation for the discretization effect for the TPL coupling constant.
This section is a preliminary report for the study on the eigenvalues of the Dirac operator. We will report the detailed studies in the independent paper in near future [68].

F.1. Perturbative analysis
Let us consider the massless staggered-Dirac operator D(x, y), where η µ (x) is the staggered phase. This eigenvalues of D are pure imaginary since it is the anti-hermitian D † (x, y) = −D(x, y): 35/41 where ψ (l) λ denotes a Dirac fermion and (l) denotes the level of the eigenvalues. We define the lowest one as l = 1. The degree of freedom of one staggered-Dirac operator is 16, and there are additional 3 color and 3 smell indices for our Dirac fermion (see Eq. (2.7)). The number of flavor (N f = 12) is realized by (4 taste's) × (3 smell's) degrees of freedom and there are 3 flavor (=smell) staggered Dirac fermions. The operator D † D is hermitian and positive definite, and it can be decomposed to the operators on the even and odd sites. In this work, we measure the eigenvalues only positive and in the even-to-even sites, and then the degeneracy of one staggered-Dirac operator (4 taste's× 4 spinor's) becomes half. The remain degrees of freedom (3 color's × 3 smell's) can be counted as the unphysical twisted momenta in the twisted directions as explained in the Sec. 4.  At the tree level, we can calculate the eigenvalue of the Dirac operator on the lattice: wherek µ denotes the momentum of the fermion field for each direction. The leading order of O(a) for low lying eigenvalues is proportional to the sum ofk µ . In our simulation there is twisted boundary condition and the non-trivial vacuum phase, and then the momentum is given by the Eq. (4.3) ask θ µ . The lowest momentum is given in the following case of (n ph µ , n ⊥ µ ) = (−1, 2) or (0, 0) for both µ = x and y, n ph µ = 0 for both µ = z and t.

36/41
The degree of degeneracy of that is 8, and the sum of the momentum is √ 18π/L. The third lowest mode is also counted, it is given by (n ph µ , n ⊥ µ ) = (−1, 2) or (0, 0) for both µ = x and µ = y, n ph µ = 0 for µ = z or t, n ph µ = −1 for µ = t or z (in same order).
The number of degrees of degeneracy of that is 8, and the sum of the momentum is √ 22π/L. Let us compare the measured value of the simulation with this tree level analysis. The total degree of the degeneracy should be multiplied by 8, since we measure the eigenvalue of the staggered fermion on even-to-even site as we explained. Figure F1 shows the first 100 eigenvalues in the β = 50, L/a = 6. There is a clear gap between l = 32 and l = 33 as we expected, although there is a large taste breaking since the lattice size is small. The ratio of the eigenvalue between first and 33rd levels is 1.328(2), and it is almost consistent with the tree level prediction 18 10 . On the other hand, we cannot see the second gap, which we expect to lie between 96th and 97th levels. The numerical value of the ratio of them is 1.110(1), and the value is completely consistent with 22 18 . Although the second gap is not clear, the eigenvalue reproduces the tree level prediction.

F.2. The eigenvalue distribution and the taste breaking
We also measure the eigenvalue for all lattice parameters in the data sets in Appendix A. We show some examples in Fig. F2, in which the error is estimated by the bootstrap analysis. The eigenvalues of high β and the large lattice size show the 4-fold degeneracy, although there is large taste breaking in the small lattice size or low β region. The data at the lowest β, β = 4.5 for L/a = 6, 12 and β = 5.5 for L/a = 20, show the inconsistency with zero, and the β dependence of the data at fixed lattice extent is smooth in whole β region. If we assume that the Banks-Casher relation, and the chiral symmetry is also preserving even in the lowest β in our simulation parameter. To see the taste breaking in our step scaling analysis in the strong coupling regime, we would like to show the eigenvalues in the continuum limit by taking the TPL renormalized coupling as a reference of the same physics: L · λ where N is the number of the fitting parameter, and in practical we choose the best fit value of N for each l and L/a. Typically, N = 4 -7 are employed in this analysis. In this analysis, the leading discretization error comes from the eigenvalue itself, which is proportional toL · λ (l) ∝ const. + O(a 2 ) if the theory lives in the deconfinement phase. The other leading contribution O(a 2 ) comes from the renormalized coupling, and the contribution comes via the tuned value of β. We take the continuum limit for 5 data points, L/a = 38/41 8, 10, 12, 16 and 20, and they can be fitted well by the linear function of (a/L) 2 in whole region. To estimate the systematic uncertainty of this continuum extrapolation we also show the quadratic extrapolation of (a/L) 2 for 6 data points included the coarsest lattice L/a = 6.  Fig. F4 The continuum extrapolation for the low lying eigenvalues for u = 2.5. Each full symbol denotes the eigenvalue at the interpolated β for each lattice size; L/a = 6, 8, 10, 12, 16 and 20 from the right to left. The shadow symbol at (a/L) 2 = 0 denotes the extrapolated value by the linear extrapolation function of (a/L) 2 for the finer 5 lattice data. The empty symbol at (a/L) 2 = −0.001 shows the extrapolated value for the quadratic function of (a/L) 2 by using all 6 data points.
As an example, we take u = 2.5, which corresponds to the region for (β, L/a) = (5.4, 6) - (6.4, 20). The taste breaking of the raw data in these region is strong. The continuum extrapolation for u = 2.5 is shown in Fig. F4. In this plot, we show the eigenvalues at continuum limit for the low lying eigenvalues; level 1 -8. The shadow symbol at (a/L) 2 = 0 denotes the extrapolated value for the linear extrapolation function of (a/L) 2 for the finer 5 lattice data. The empty symbol at (a/L) 2 = −0.001 shows the extrapolated value for the quadratic function of (a/L) 2 by using all 6 data points. The difference between these two kinds of the extrapolation can be identified as the systematic uncertainty. Including the systematic error, we find that the breaking of the level 1 -4 becomes mild at the continuum limit even in the strong coupling regime. We also consider the order a 4 effect in the running coupling study (Fig. 18). The recoveries of the taste breaking in the continuum limit using the extrapolation function at the same order would show an indirect evidence that we have considered sufficient order in a 2 in taking the continuum extrapolation of the TPL coupling constant.