A novel scheme for the wave function renormalization of the composite operators

We propose a novel renormalization scheme for the hadronic operators. The renormalization factor of the operator in this scheme is normalized by the correlation function at tree level in coordinate space. If we focus on the pseudo scalar operator, then its renormalization factor is related to the mass renormalization factor of the fermion through the partially conserved axial-vector current (PCAC) relation. Using the renormalization factor for the pseudo scalar operator in our scheme, we obtain the mass anomalous dimension of the SU(3) gauge theory coupled to N_f=12 massless fundamental fermions, which has an infrared fixed point (IRFP). The mass anomalous dimension at the IRFP is estimated as gamma_m^*= 0.044_{-0.024}^{+0.025} (stat.)_{-0.032}^{+0.057} (syst.).

We propose a novel renormalization scheme for the hadronic operators. The renormalization factor of the operator in this scheme is normalized by the correlation function at tree level. If we focus on the pseudo scalar operator, then its renormalization factor is related to the mass renormalization factor of the fermion through the partially conserved axial-vector current (PCAC) relation. Using the renormalization factor for the pseudo scalar operator in our scheme, we obtain the mass anomalous dimension of the SU(3) gauge theory coupled to N f = 12 massless fundamental fermions, which has an infrared fixed point (IRFP). The mass anomalous dimension at the IRFP is estimated as γ * m = 0.044 +0.025

Introduction
Lattice gauge theory provides a regularization method for the gauge theory. To regulate the theory, we introduce a lattice spacing (a) as a ultraviolet (UV) cutoff and a finite lattice extent as an infrared (IR) cutoff. For the lattice gauge theory, there are several useful renormalization schemes for the gauge coupling constant, e.g. the Schrödinger functional (SF) scheme [1], the potential scheme [2], the Wilson loop scheme [3], the twisted Polyakov loop (TPL) scheme [4,5], the Wilson flow (Yang-Mills gradient flow) scheme [6,7] etc. A variety of renormalization scheme for composite operators has also been given, e.g. the SF scheme [8] and RI-MOM scheme [9] and so on. Concerning the fermion mass renormalization, the quark mass renormalization factor is related to that of the pseudo scalar operator because of the partially conserved axial-vector current (PCAC) relation.
In this paper, we propose a novel scheme for the composite operators. The basic idea is to normalize the renormalization factor using the tree level correlation function of the operator. One of the practical advantages of this scheme is that, in principle, a special lattice setup (such as a special boundary condition etc.) is not necessary.
We also apply this renormalization scheme to derive the anomalous dimension of the pseudo scalar operator for the SU(3) gauge theory coupled to N f = 12 massless fermions. In our previous work [5], we investigated the running coupling constant of this theory using the twisted Polyakov loop scheme from the perturbative region to the IR region. We found the growth of the renormalized coupling halts in the IR region, which verifies that the infrared fixed point (IRFP) exists in this theory. At the IRFP, we expect that an interactive conformal field theory is realized. Note that the lattice gauge action is defined at the Gaussian UV fixed point and we do not know the explicit form of the action of such interactive IR conformal theory. However, we expect that the theory is sufficiently close to the IRFP in the region, where the coupling constant does not show the growth when the energy scale changes. Here we declare that the theory on the lattice realizes the conformal fixed point theory.
Conformal fixed points are the most important object in the quantum field theories. At the conformal fixed point the critical exponents (e.g. the anomalous dimension of the operators) are the scheme independent quantities, and these exponents classify the universality class of the quantum field theories. Among the several critical exponents, the one that is related to the relevant operator is crucial to define the IR field theory. In this paper, we determine numerically the universal mass anomalous dimension of the interactive conformal field theory.
Recently, several methods to obtain the mass anomalous dimension for the conformal gauge theory realized at the IRFP in many-flavor SU(N c ) gauge theories have been proposed. The step scaling is one of the methods based on the renormalization group for the finite scaling, and this method can be applied to the non-conformal field theories [1]. The other method is to use the hyperscaling [10] for the mass deformed conformal gauge theory. The application of the hyperscaling on the lattice was pointed out by M. A. Luty and F. Sannino [11,12], and the detailed practical discussions of mass deformed conformal gauge theory have been shown by L. Del Debbio and R. Zwicky [13]. This method is based on the assumption of the existence of the interactive conformal field theory, where the scaling of the operators is different from the Gaussian (canonical) case. The mass anomalous dimension is derived by the fit of the mass spectrum of hadronic state or the chiral condensate in a small mass region. A similar method to estimate the mass anomalous dimension using the fit for the massless SU(3) gauge theory has also been proposed [14,15]. They utilized the massless fermion, and from the scaling of (1/L), where L is a finite lattice extent, they estimated the mass anomalous dimension. The independent method has been suggested in the paper [16]. They assumed that the correlation function in the finite volume around the IRFP became the Yukawa-type function and derived the anomalous dimension from the fit.
In this work, we obtain the mass anomalous dimension at the IRFP using the step scaling method. In our numerical simulation, we introduce the twisted boundary condition for both gauge field and the fermion field. The boundary condition kills the zero mode contribution and regularizes the fermion matrix even in the massless case. Thus we carry out the simulation using exactly massless fermions. Several independent groups have been investigating the mass anomalous dimension of the SU(3) gauge theory coupled to N f = 12 fermions [15,17] - [19]. The works [17] - [19] are based on the hyperscaling method for mass deformed conformal gauge theory applied to the simulation with massive fermions. In the paper [15], they utilize the (approximately) massless fermion, and derive the universal mass anomalous dimension in the infinite volume limit using the hyperscaling for the Dirac eigenmodes. This work is the first study on the mass anomalous dimension for the SU(3) N f = 12 massless gauge theory using the step scaling method. We expect that the value of the mass anomalous dimension at the IRFP is independent of the derivation. This paper is organized as follows: In Sec. 2, we give the definition of a novel renormalization scheme for the composite operators. The renormalization factor of the pseudo scalar operator is related to the fermion mass renormalization factor, therefore in the rest of the paper we focus on the pseudo scalar operator and determine the fermion mass anomalous dimension. In Sec. 3, we show the strategy to obtain the mass anomalous dimension using the step scaling method and give a definition of the step scaling function in our scheme.
Note that there are two definitions of the step scaling function in our scheme, since there are two scales in the observable for the renormalization condition. In Sec. 4, we explain our numerical simulation setup. We compute the correlation function at tree level, which is needed to define the renormalization factor in our new scheme in Sec. 5 and Appendix A. We determine the mass anomalous dimension at the IRFP of the SU(3) gauge theory coupled to N f = 12 massless fermions in Sec. 6. We find that the mass anomalous dimension at the IRFP is given by where the first systematic error comes from the uncertainty of the continuum extrapolation while the second one comes from the uncertainty of the value of the coupling constant at the IRFP. We discuss the comparison with the other works in Sec. 7. We conduct a discussion about the promising methods to studying such universal quantity around the IRFP using the lattice simulation in Sec. 8.

A novel renormalization scheme for the anomalous dimension
We give a new renormalization scheme of some arbitrary composite operator (H). In renormalizable theories, a nonperturbative renormalized coupling constant can be defined by amplitudes of the observables. The SU(3) gauge theory coupled to a small number (N f ≤ 16) of fundamental fermions is asymptotically free and it is described by two kinds of parameter: the gauge coupling and the mass parameter of the fermions. The renormalization factor can be defined by the correlator of the bare operator (H), To obtain the finite renormalized value of the correlator, we introduce a nonperturbative renormalization factor (Z H ) as follow: Here C R H denotes a renormalized correlation function and it is finite. On the other hand, the renormalization factor Z H and the nonperturbative bare correlation function diverge, and on the right hand side these divergences are canceled each other.
We introduce the renormalization condition on the renormalized correlator, in which the renormalized correlator is equal to the tree level amplitude: The renormaliation factor of the composite operator is thus defined by at some fixed propagation length (t). Thus the factor Z H is normalized by the tree level value for each propagation length. On the lattice, the nonperturbative bare correlation function (C H (t)) is calculated by lattice numerical simulations. The correlation function on the lattice depends on the propagation time (t/a), the bare coupling constant (g 0 ), the bare mass (m 0 ) and the lattice 3/26 size(L/a, T /a). Thus it is denoted by C H (g 0 , m 0 a, t/a, T /a, L/a) on the lattice. We fix the ratio between the temporal and the spatial lattice extents, and identify an inverse of lattice spatial extent (1/L) as a renormalization scale (µ). Then, we give the definition of Z factor in this scheme on the lattice, Z H (g 0 , m 0 a, t/a, a/L) = C tree H (g 0 , m 0 a, t/a, a/L) C H (g 0 , m 0 a, t/a, a/L) . (2.5) Let us define these nonperturbative renormalized parameters at the energy scale (µ) as g 2 (µ) ≡ Z g g 2 0 andm(µ) ≡ Z m m 0 . Here Z g and Z m denote nonperturbative renormalization factors for each parameter. These factors are functions of dimensionless parameters g 0 ,m 0 a and aµ, so that they can be written by Z g = Z g (g 0 , m 0 a, aµ) and Z m = Z m (g 0 , m 0 a, aµ) 1 . In QCD, the renormalized mass is also defined through the partially conserved axial-vector current (PCAC) relation: where A R and P R denote the renormalized axial-vector current (A Rµ (x) = Z Aψ (x)γ µ γ 5 ψ(x)) and pseudo scalar operator (P R (x) = Z Pψ (x)γ 5 ψ) respectively. Here we introduce the renormalization factors for these operators. Note that Z A (g 0 , m 0 ) is scale independent because the axial current is renormalized through current algebra. Thus the PCAC relation gives the relationship between the mass renormalization factor (Z m ) and the renormalization factor (Z P ) of the pseudo scalar operator as, at each renormalization scale. The anomalous dimension of the dimensionless running mass of the fermions, µdm/dµ = γ m (ḡ,m)m, (2.8) can be calculated from the scale dependence of the pseudo scalar renormalization factor Z P as In this paper, we study on the massless fermion theory, and measure the renormalization factor of the pseudo scalar operator. One of the practical advantages of our renormalization scheme is that a special lattice setup, such as a special boundary condition, is not needed.

Step scaling function
Let us consider the scale dependence of the renormalization factor Z P (Eq. (2.9)). First, we introduce the discrete mass step scaling function from the factor Z P .
In our renormalization scheme, the renormalization factor Z P has two independent scales, the propagation time (t) and lattice temporal size (T ) 2 . To see the scale dependence of the factor Z P , there are two definitions of the scaling function. One definition of the scaling function is given by the factor Z P at a fixed ratio r = t/T , which we call as "fixed r" definition. Here r takes a value 0 < r ≤ 1/2 because of the periodic boundary condition on the lattice. We obtain the scale dependence of the factor Z P when we change both the physical propagation length and lattice size together. The other definition is given by the factor Z P at a fixed t. We change the dimensionless quantity r and the temporal physical lattice extent T in the latter definition. In the former definition, the renormalization scale is parametrized by the lattice spatial size. In the latter case, if T ≪ L and t > L, then the renormalization scale is given by 1/t. In this paper, we use the former "fixed r" definition. Now, the factor Z P in Eq. (2.5) depends only on the bare coupling constant and the lattice spatial size: Z P (g 0 , a/L) in a fixed r scheme. From now on we use β to denote the bare coupling constant and β = 6/g 2 0 . We give a comment on a choice of the parameter r in "fixed r" scheme. If we choose small r, the lattice data might suffer from a large discretization effect. On the other hand, in the case of large r, the signal of correlators of the hadronic operators might become noisy except for the lightest state. Practically, we have to search for the optimal range of r. Note that at the fixed point, the anomalous dimension is a scheme independent quantity, so that it should be independent of r. We discuss r dependence of our result in Sec. 6.3. Now, we give a brief review of the strategy to obtain the mass anomalous dimension using the step scaling method. The idea is established by ALPHA collaboration [8]. Figure 1 shows a schematic picture which describes the strategy of step scaling for the renormalized coupling constant and the renormalization factor of operators. The top panel shows the step scaling Fig. 1 The strategies of the step scaling method for the coupling constant and the renormalization factor for the operator P . We measure the growth ratio of each quantity at fixed β.

5/26
for the renormalized coupling. To obtain the scale dependence of the renormalized coupling in a renormalization scheme, we measure the growth ratio of renormalized coupling when the lattice extent becomes s times with fixed value of bare coupling constant. Practically we carry out the following procedures. First, we choose a value of renormalized coupling constant u = g 2 R (1/L) and tune the value of β to realize u for each lattice size. Next, we measure the renormalized coupling constant with the tuned value of β at the larger lattice sL/a. The renormalized coupling constant on the larger lattice is called the discrete step scaling function: Σ(u). Here s is the step scaling parameter (1 < s). Finally we take a continuum limit of the discrete step scaling function: σ(u) = lim a→0 Σ(u) = g 2 R (µ = 1/sL). The growth ratio of the renormalized coupling (σ(u)/u) essentially gives a discrete beta function.
To obtain the scale dependence of the renormalization factor of a operator, we measure the growth ratio of the factor Z P . If the operator is a pseudo scalar operator, it is called the discrete mass step scaling function because of the relationship Eq.(2.7). The explicit definition of the mass step scaling function is given by The mass step scaling function on the lattice includes the discretization error. To remove it, we take the continuum limit (a → 0) keeping the renormalized coupling (u = g 2 R (1/L)) constant. (3.2) In the continuum limit, this mass step scaling function is related to the mass anomalous dimension as, where γ m (x) and β(x) denote the mass anomalous dimension and the beta function respectively, and d 0 and b 0 denote coefficients of them in 1-loop order. This relation becomes simple when the theory is conformal: and we can estimate the anomalous dimension at the fixed point with the following equation: Here u * denotes the fixed point coupling constant. Note that in Eq. (3.2), there is a freedom of the choice of the renormalization scheme concerning the input renormalized coupling constant. We use the set of the bare coupling constant and the lattice size to realize the input renormalized coupling constant (u) with a renormalization scheme for the gauge coupling constant. The energy scale is defined by the input renormalized coupling, and the energy dependence of the mass step scaling function comes through the renormalized coupling constant. We can use any combinations of the renormalization schemes for the renormalized coupling and the wave function renormalization. Generally, the value of σ P (u) and the mass anomalous dimension depends on 6/26 the choice of the renormalization schemes. At the fixed point, although the value of u * depends on the renormalization scheme, the mass anomalous dimension is independent of the renormalization schemes of both the mass and the coupling constant.

Simulation setup
The gauge configurations are generated by the Hybrid Monte Carlo algorithm, and we use the Wilson gauge and the naive staggered fermion actions. We introduce the twisted boundary conditions for x, y directions and impose the usual periodic boundary condition for z, t directions, which is the same setup with our previous work [5]. Because of the twisted boundary conditions the fermion determinant is regularized even in the massless case, so that we carry out an exact massless simulation to generate these configurations. The simulations are carried out with several lattice sizes (L/a = 6, 8, 10, 12, 16 and 20) at the fixed point of the renormalized gauge coupling in the TPL scheme [5]. In this simulation, we fix the ratio of temporal and spatial directions: T /a = 2L/a 3 . We use the tuned value of β where the TPL coupling is the fixed point value for each (L/a) 4 lattices. We neglect the possibility of induced scale violation coming from the change of lattice volume (L/a) 4 → 2(L/a) 4 since we carefully take the continuum limit. We measure the pseudo scalar correlator for 30, 000-80, 000 trajectories for each (β, L/a) combination. We estimate the statistical error using bootstrap method. Now, we explain the detailed values of β in our simulation. In this paper, we focus on the mass anomalous dimension at the IRFP. In our previous paper [5], we found that the existence of the IRFP at g * 2 TPL = 2.69 ± 0.14(stat.) +0 −0.16 (sys.), (4.1) in the TPL scheme. The parameter sets which realize the fixed point coupling on the lattice is shown in Table. 1. We generate the configurations using these parameters on (L/a) 3 Table 2 Additional simulation parameters.
The pseudo scalar correlator can be presented by the fermion propagators (S(t, x)) where γ 5 ⊗ γ 5 specifies the spin and flavor structure. We measure the pseudo scalar correlator using the point source at t = 0. We construct the Dirac field by the staggered fermion(χ) on the hypercubic space-time. The pseudo scalar correlator calculated by the two point function has the value at even temporal sites: Practically, we introduce a tiny bare fermion mass ma = 10 −5 -10 −6 for the measurement of the correlators. To check the smallness of this bare quark mass rather than the twisted momentum even in the strong coupling region, we have changed the mass to ma = 10 −7 and confirmed that the effect of the mass is neglegible.

Calculus of the correlator at tree level
We compute the correlator at the tree level to normalize the renormalization factor in the renormalization condition Eq. (2.3). The correlation function of the pseudo scalar at the tree level corresponds to the correlation function of the two free fermions. We can calculate it using the vacuum configurations. There are three possible choices of the vacuum configurations in the case of SU(3) gauge theory, since the pure SU(N c ) gauge theory has Z Nc global symmetric degenerated classical vacua at U µ = exp(2πiθ µ /N c ), where θ µ = 0, 1, · · · , N c − 1 for each direction. According to the semi-classical analysis of the 1− loop effective potential (See Sec. 4 in Ref. [5]), the 3 4 −fold degenerated vacua, where the Polyakov loop in z, t directions has a nontrivial complex phase (exp(±2πi/3)) are chosen in our lattice setup. Therefore we use that the vacuum configuraions, U µ = exp(±2πi/3T )I for µ = z, t, to derive the correlator at the tree level. For x and y directions, the effective potential does not depend on the choice of the θ µ . We use the simple constant configuration U µ = I. The data of the pseudo scalar correlator for each lattice size using this vacuum configuration is shown in Table B1 in Appendix B. The data can be fitted by cosh[ω(t −T /2)] well in longt region (See Appendix A), where the hatted symbol denotes the quantities in the lattice units. Let us consider the meaning of "ω" in the massless fermion case. Figure 2 shows the values of the fitted parameter "ω" for several lattice extents (L/a) 3 The value of ω when we solve the function C(t) = c 0 cosh(ω(t −T /2)) using two independent data aroundt =T /2. The dot line denotes ω 0 = 2E 0 , where E 0 denotes the lowest energy of single free fermion given by Eq.(5.3).
and 40. For each T /L, the parameter "ω", which we call the effective mass, is proportional to a/L. At the long distance, only the lowest energy mode must survive, and we expect that the correlation function can be approximated as C(t) ∼ e −2E0t , where E 0 is the lowest energy of single fermion. It is realized at the lowest energy state those four-dimensional momentum is zero: Thus the lowest energy is obtained by the sum of lowest spatial discrete momenta. In our lattice setup the momentum for each direction is given in Eq. (22) in the paper [5]:p where n µ = 0, 1, · · · ,L/2 − 1 and m µ = 0, 1, · · · , N c − 1 with (m ⊥ x , m ⊥ y ) = (0, 0). The lowest energy of single fermion is analytically calculated as the following, The dot line in Fig. 2 denotes the line of ω 0 = 2E 0 . We find that around t/a = T /2a with T /L ≥ 4 only the lowest mode remains. In this paper, we uese T /L = 2 lattices, so that there are still some contributions from the second lowest energy and the higher modes in the correlator at the tree level. However, the lattice data in Fig. 2 shows that the effective mass ω vanishes in the continuum limit even if such modes remains.

Mass anomalous dimension at the IRFP
We measure the mass anomalous dimension at the IRFP for the SU(3) N f = 12 gauge theory using the step scaling method. When we take the continuum limit of the mass step scaling 9/26 function, we use the TPL coupling constant at the IRFP as an input renormalized coupling. In the paper [5], the value of TPL coupling at the IRFP is determined We take the mass anomalous dimension at the central value, (g * 2 TPL = 2.686) as a central analysis using the s = 1.5 step scaling for the "fixed r" scheme with r = 1/2. We estimate the systematic error by taking the discrepancy among several kinds of the continuum extrapolation. We also derive γ * m at the lower bound of the fixed point coupling (g * 2 TPL = 2.475) and the upper bound of the coupling (g * 2 TPL = 2.823) to include the systematic uncertainty coming from the value of the fixed point coupling. We also show the results of the s = 2 step scaling and the dependence on the scheme parameter r.
6.1. Result of s = 1.5 step scaling function in r = 1/2 scheme First, we compute the mass anomalous dimension at the IRFP using the s = 1.5 step scaling. We show the renormalization factor Z P (β, a/L, t/a) in our scheme in Fig. 3  raw data of the factor Z P (t/a) for each lattice setup in Tables B2 -B7 in Appendix B. We found that each Z P (t/a) has a different slope between the long t/a and short t/a regions. That means that the contributing effective mass depends on the distance. The data in the short propagation length (t/a ≤ 6) for each lattice size seem to have a large discretization effect. To reduce the effect, we choose r = t/T = 1/2 scheme in our central analysis.
To carry out the s = 1.5 step scaling, the interpolations of the data in β and L/a are necessary. Figure 4 shows the β dependence of Z P (β, L/a) at t/a = T /2a. We find that the 10  Finally, Fig. 5 shows the mass step scaling function Σ(β, a/L; s = 1.5) on the lattice with the scheme parameter r = 1/2. We find that the data L/a = 6 data suffers from a larger discretization error. The discretization effects arise from two sources. One is a discretization effect of the renormalized coupling due to the tuned value of β. The other comes from a discretization effect of the pseudo scalar correlator. As we shown in Fig. 3, since the data in the short propagation range (t/a ≤ 6) seems to have a large discretization, the latter is the dominant source of the large scaling violation.
Since the fit including the data at L/a = 6 has a large chi-square, we drop the data at L/a = 6 from the continuum extrapolation 4 . Figure 6 shows the finer three lattice data and The error denotes the statistical one, which is estimated by the bootstrap method. We also carry out two different kinds of extrapolation. One is the three-point constant extrapolation (the blue dashed line in Fig. 6) and the other is the two-point linear extrapolation (the violet dashed line in Fig. 6). The smallest value of γ * m is given by the three-point constant extrapolation for u = 2.475, and the largest one is given by the two-point linear extrapolation for u = 2.686. Each value of γ * m is 0.013 and 0.102 respectively. We estimate 4 The same situation happened in the previous work for the running coupling constant [5]. At that time, since L/a = 4 data of g 2 TPL on the lattice suffer from large discretization effects, we dropped the data from the continuum extrapolation. where the first systematic error comes from the uncertainty of the continuum extrapolation while the second one comes from u * dependence. Note that the corresponding degrees of freedom for each three-point constant, three-point linear and two-point linear extrapolations are 2, 1 and 0 respectively. The extrapolation with the small degree of freedom might srtongly suffer from the statistical fluctuation. Furthermore, there is a signal that the finer lattice data would give a large discrepancy from the unity line, thus it might give a large anomalous dimension. The further study including the larger lattices is necessary to give a conclusive result.

Step scaling parameter (s) dependence
We also show the result of the step scaling with s = 2. The advantage of the s = 2 step scaling is that we do not need interpolations in β and L/a, and the signal of the growth ratio of the factor Z becomes clear. On the other hand, the disadvantage of this is that we need the large lattice setup with the fixed distance from the continuum limit. These differences are just technical matters, and we can make the consistency check when we carry out the step scaling with several values of s. Note that the anomalous dimension within the same renormalization scheme is independent of s, although the mass step scaling function depends on the choice of s.   Figure 7 shows the mass step scaling function in the case of s = 2 in r = 1/2 scheme. Again we find that there is a large discretization error in L/a = 6, therefore we drop the data. The remained data points are only two, so that we take the average of these two points for each input u. The degree of freedom of the continuum extrapolation is the same with the 13/26 central analysis in the s = 1.5 step scaling analysis. The anomalous dimension is obtained γ * m = 0.028 ± 0.006 for u = 2.475, γ * m = 0.023 ± 0.007 for u = 2.686 and γ * m = 0.034 +0.007 −0.008 for u = 2.823. These results are consistent with the ones of s = 1.5 analysis within 1σ. That is an indirect check for the effects of the β and L/a interpolations in the s = 1.5 step scaling.

Scheme parameter (r) dependence
We also show the result in the scheme with r = 1/3. Changing r corresponds to the change of renormalization condition, so that it gives a different renormalization scheme with the same lattice setup. If the theory is not at the fixed point, the value of the anomalous dimension depends on the renormalization scheme, namely the choice of r. However at the fixed point, it should be independent of the choice of the renormalization scheme.
In r = 1/3 scheme how to estimate the value of correlator at non-integer t/a in several lattice sizes is a problem. As discussed in Appendix A, we expect that the correlation function in the finite box can be described by the exponential functions even though the theory is conformal. Since we impose the periodic boundary condition in the temporal direction, the correlator is proportional to the linear combination of the cosh function of several energy Here ω i is the effective mass for each mode. We assume that our data can be fitted by a single cosh function in the small range oft. Figure 8 shows the examples of thet interpolation fort =T /3. We take the following fit range for each lattice size: (4 ≤ t/a ≤ 6) and (10 ≤ t/a ≤ 12), (6 ≤ t/a ≤ 8) and (12 ≤ t/a ≤ 14) , (10 ≤ t/a ≤ 12) and (20 ≤ t/a ≤ 22) and (12 ≤ t/a ≤ 14) and (26 ≤ t/a ≤ 28) for L/a = 8, 10, 16 and 20 respectively. The blue curve in Fig. 8 denotes the fit function. Note that essentially the number of independent data points is two, so that we solve the equation C P S (t) = a cosh(b(t −T /2)) to determine the fit parameters (a and b).
We also carry out the same interpolation for the tree level correlators, and repeat the same analysis as we shown in Sec. 6.1. Figure 9 shows the mass step scaling function for each input value of u. We find that there are the large scale violations in particular L/a = 6 and 8 data. As we explained, the origin of these scale violations comes from two sources: the discretization error of the input renormalized coupling constant and the one of the correlation function. Both Fig. 5 and Fig. 9 have the same discretization error coming from the scale violation of the coupling constant. The difference of Σ(β, a/L) between these two plots for each data point shows the discretization error of the pseudo scalar correlator. The data at L/a = 6 and L/a = 8 include the data of Z P (t/a) at t/a ≤ 6 in Fig. 3. We expect that such data in the short propagation length gives the large scaling violation.
We estimate the value of step scaling function in the continuum limit using the constant extrapolation for two finer lattice data. Thus the degree of freedom of the continuum extrapolation is the same with the central analysis in r = 1/2 scheme again. The anomalous dimension is given by γ * m = 0.020 ± 0.007 for u = 2.475, γ * m = 0.037 ± 0.008 for u = 2.686 and γ * m = 0.050 ± 0.008 for u = 2.823. These results are also consistent with the result (6.3) within 1σ. It is an evidence that the anomalous dimension at the IRFP shows the universal property. 14/26

Discussion
Our result of the mass anomalous dimension at the IRFP is where the first systematic error comes from the uncertainty of the continuum extrapolation while the second one comes from u * dependence. Let us compare our result with the other predictions. Figure 10 shows the values of the mass anomalous dimension in other literatures. The perturbative results [20][21][22] [17], [18], [19], [15] and our result. Note that in the papers [18,19] there is no " * " on the gamma in their own papers.
the perturbative prediction in the case of N f = 12, but it is the same order with N f = 16 case. Furthermore, it is also important to compare our result with the other nonperturbative lattice studies. There are four papers to estimate the mass anomalous dimension in this theory using the mass deformation method. The most reliable result, γ * m = 0.32 (3), is given by Cheng et. al. in the paper [15]. The hyperscaling (volume-scaling) of the Dirac eigenmodes for several values of β has been investigated in the (approximately) massless limit. The results show a reliable scaling behavior and they can estimate both the fixed point value of β in 16/26 their lattice gauge action and the universal mass anomalous dimension in IR limit. The discrepancy between our result and their result might come from an insufficient estimation of the systematic uncertainty coming from the continuum extrapolation in our analysis.
The papers [17] and [18] use a part of the data in the paper [23], and fit these data using the hyperscaling, which is the Miransky scaling [10] on the lattice. Orginally, the paper [23] shows that the data can be fitted by the weakly chiral symmetry broken hypothesis better than the conformal hypothesis. However, the paper [17] shows if we fit only the largest lattice size data to avoid serious finite volume effects, the conformal hypothesis also works well. The universal fit using the hyperscaling for several hadronic spectrum gives the mass anomalous dimension: γ * m = 0.403 (13). The paper [18] shows the finite-size scaling using the same data, and the anomalous dimension of the pseudo scalar operator from the mass spectrum is given by γ m = 0.35 (23). In the paper [19] by LatKMI collaboration, they also uses the same method, and fit their own data. The anomalous dimension in paper [19] is given as γ m ∼ 0.4-0.5. Here, we should mention that in the papers [18,19] the symbol " * " is not added to the anomalous dimension, which denotes a symbol of the quantities at the IRFP. Actually, LatKMI also has been studying the mass anomalous dimension in the case of N f = 16 fermions using the same strategy [24]. However, the result of the paper [24] using the hyperscaling also shows a scaling behavior within the wide range of the larger mass anomalous dimension, and do not show the perturbative converged value of the mass anomalous dimension is γ * m ∼ 0.026. If it is not the anomalous dimension at the fixed point, then it depends on the renormalization scheme and the discrepancy is not a problem.
In fact, the SU(2) gauge theory coupled to N f = 2 adjoint fermion is also known as an IR conformal field theory. Several independent collaborations have been deriving the mass anomalous dimension at the IRFP. The step scaling method using the SF scheme gives the predictions 0.05 ≤ γ * m ≤ 0.56 in paper [25] and γ * m = 0.31 (6) in paper [26] respectively. The hyperscaling for the string tension, Meson spectrum and the mode number of the Dirac operator give γ * m = 0.22(6) (Ref. [27]), 0.05 ≤ γ * m ≤ 0.20 (Ref. [28]) and γ * m = 0.51(16) (Ref. [29]), and γ * m = 0.371(20) (Ref. [30]). There are somewhat consistent with each other, while some values have a large errorbar and there are also unknown systematic errors. The paper [25] determined the critical value of β around β ∼ 2.25 using the step scaling method in the SF scheme. Also in the paper [31], the Creutz ratio does not run around β = 2.25 and that is a signal of the fixed point of the Wilson loop coupling. 5 . Then in the papers [27] - [29], they derive the mass anomalous dimension using the hyperscaling with tuned value of β = 2.25. We consider that such tuned value of β to realize the IRFP is necessary to obtain the universal anomalous dimension using the hyperscaling for the mass deformed gauge theory.

Summary
We propose a new renormalization scheme for the composite operators. In this renormalization scheme, the correlator of the pseudo scalar operator satisfies the "tree level renormalization condition", in which the renormalized value is equal to the tree level amplitude. In this scheme, the different propagation length corresponds to the different renormalization schemes, and we choose the suitable length in practical simulations.
Furthermore we study the mass step scaling function for the SU(3) N f = 12 massless gauge theory using this renormalization scheme. Using the PCAC relation, the mass renormalization factor is related to the renormalization factor of the pseudo scalar operator. We actually measure the renormalization factor of the pseudo scalar operator, and directly derive the mass anomalous dimension. This work is the first study of the derivation of the mass anomalous dimension at the IRFP using the step scaling method.
Our result of the mass anomalous dimension at the IRFP of this theory with the scheme parameter r = 1/2 and the step scaling size s = 1.5 is where the first systematic error comes from the uncertainty of the continuum extrapolation while the second one comes from u * dependence. We also investigate the step scaling parameter (s) dependence and the scheme parameter (r) dependence of γ * m . The results with the different choices of s and r are consistent with each other within 1σ discrepancy. Note that in the current analysis the continuum extrapolation is done with small degrees of freedom, and the result might be strongly affected by the statistical fluctuation. Furthermore, there is a signal that the finer lattice data would give a large discrepancy from the unity line, thus it might give a large anomalous dimension. Further careful estimation of the systematic uncertainty from the continuum extrapolation might be important. We will report the conclusive results including the larger lattice ((L/a) 3 × (T /a) = 24 3 × 48) simulation in the forthcoming paper.
If there is no other relevant operator, then the renormalization group flows of the SU(3) N f = 12 gauge theory are governed by the two dimensional theory spaces whose coordinates are the fermion mass and the gauge coupling constant (See: Fig. 11). The universal quantities γ γ Fig. 11 The theory space for the SU(3) N f = 12 gauge theory.
to characterize the IRFP are the critical exponent of the beta function (γ * g ) and the mass anomalous dimension (γ * m ). We have investigated the renormalization group flow on the massless line. We also derived in our previous work [5]. We determine γ * g and γ * m using the TPL scheme for renormalized coupling constant and the new scheme for the fermion mass respectively. Changing the renormalization schemes corresponds to the coordinate transformation of the theory spaces. The existence of IRFP is independent object to the coordinate transformation. The values of γ * g and γ * m are also universal, since they are the eigenvalues of two linearized β functions around the IRFP.
We compare our result with other lattice studies. All other studies has been done based on the scaling law. There is a large difference between our result and their, but some works utilize the mass deformed theory without the tuning of β. We consider that an insufficient parameter tuning for the hyperscaling could be a reason of the discrepancy. The determination of γ * m using the hyperscaling method works well only if the renormalization group flow reach as the vicinity of the IRFP such as the solid curve in Fig. 11. Since we do not know the action at the IRFP, we introduce the fermion mass term in the lattice gauge action at the UV gaussian fixed point. Around the gaussian fixed point, the mass term is relevant operator and the gauge coupling is marginal, so that generally the renormalization group flow goes away from the massless axis such as the dotted curve in Fig. 11. If it happens, the renormalization group flow reaches the renormalized trajectory (RT in Fig. 11) where it is far away from the IRFP. The anomalous dimension changes along the renormalization group flow even on the RT. If there is no scale invariance and the mass anomalous dimension is not the one at the IRFP, then generally the value of γ m depends on the renormalization scheme.
To find both the IRFP and to obtain the universal mass anomalous dimension needs two independent observables. It is impossible to do both only using the hyperscaling for the mass deformed gauge theory in two parameter spaces (β, m).
One of the most promising methods to find the IRFP is the step scaling method for the renormalized coupling constant. First we have to find the IRFP of renormalization group flow using the method in a renormalization scheme. Next derive the mass anomalous dimension at the IRFP using the tuned lattice setup to realize the IRFP using the step scaling or the hyperscaling or the other method. The hyperscaling for the mass deformed conformal gauge theory is a powerful method to obtain the precise value of the anomalous dimension. However, we would like to emphasize the tuning of the lattice parameters (β, L/a and the fermion bare mass) in order to realize the vicinity of the IRFP is important to obtain the universal quantity.
On the other hand, the paper [15] shows the universal mass anomalous dimension from the Dirac eigenmodes using the (approximately) massless fermions with several values of β and the lattice sizes. In the massless limit, since we expect that there is only one parameter (g) in the theory space, the universal behavior appears by measuring one observable. To obtain only γ * m , the method looks promising. The discrepancy between our result and their result might come from an insufficient estimation of the systematic uncertainty coming from the continuum extrapolation in our analysis.
The future direction within our method is the simulation with the larger lattice size ((L/a) 3 × (T /a) = 24 3 × 48) to give a conclusive value of the critical exponent at the IRFP. Actually, in the present analysis the degree of freedom of the continuum extrapolation is only one, so that we did not estimate the systematic uncertainty coming from this procedure. We will report the larger lattice data in forthcoming paper. To measure the wave 19/26 function renormalization factor for other hadronic operators and to investigate the universal scaling behavior are also interesting. Furthermore, if we take the continuum limit carefully, then a study with the different lattice setup using the tuned values of β to realize g * 2 TPL (Table 1) would be promising to derive the anomalous dimension at the IRFP. That must be a nontrivial check for the universality using the lattice simulations 6 .

Acknowledgements
The idea of this novel scheme was suggested to us by T. Onogi. We would like to thank him and H. Ikeda for useful discussions and Y. Taniguchi for providing us with notes on the free fermion correlator in the infinite volume in Appendix A. The simulation codes was originally developed by H. Matsufuru in our previous works, we also thank him. We thank S. Hashimoto, Y. Iwasaki, M. Lüscher, A. Patella, F. Sannino, H. Suzuki, N. Yamada and S. Yamaguchi for giving useful comments and discussions. We also would like to thank A. Irie for making Figs. 1 and 11. Numerical simulation was carried out on Hitachi SR16000 at YITP, Kyoto University, NEC SX-8R at RCNP, Osaka University, and Hitachi SR16000 and IBM System Blue Gene Solution at KEK under its Large-Scale Simulation Program (No. 12/13-16), as well as on the GPU cluster at Osaka University. We acknowledge Japan Lattice Data Grid for data transfer and storage. E.I. is supported in part by Strategic Programs for Innovative Research (SPIRE) Field 5.

A. The shape of the correlation function
We give a comment on the functional form of the pseudo scalar correlator. If the theory is conformal, then in the infinite volume the correlation function of an operator O shows the power function, where ∆ O denotes the conformal dimension of the operator, which is given by the sum of the canonical and anomalous dimensions of the operator respectively. However, it is hard to fit the lattice data using the power function. We consider that there are two reasons why the data in our simulation shows cosh behavior not power law. First, we consider that the free massive fermion theory. In the continuum limit with infinite volume, we can calculate the correlation function of the pseudo scalar operator in this theory. Let us consider the correlation function of the pseudo scalar operator (P ). The correlation function can be obtained using two free fermion propagators where p, k denote four-dimensional momenta andp denotes a magnitude of three-dimensional momentum. In the case of the massive fermion, it is written by where K 2 (2mt) is the modified Bessel function. In the case of the massless fermion, the theory is a conformal free theory. The correlation function is given by as expected by the dimensional analysis.  Figure A1 shows the shape of the G(t, m) with the mass m = 0.80, 0.10, 0.01. If the mass is large, then the correlation function G(t, m) can be described by e −2mt . On the other hand, if the mass is a quite small, the G(t, m) reproduces the massless power correlation function. In the middle range of the mass (m ∼ 0.1), the correlation function can be described by the power function (Eq. (A4)) in the short t range, while it is consistent with the exponential function e −2mt in the long t region. If the mass becomes small, the available range of the power function fit becomes broad. That is related to the fact that the massless free fermion theory is an UV fixed point, so that a short range behavior describes a conformal behavior.
On the lattice, since we introduce other two scales: the finite lattice size and a lattice spacing, the discussion becomes complicated. In our simulation, the twisted boundary condition shifts the zero momentum to nonzero values. The shifted momentum plays a similar role to the mass.  Figure A2 shows the data of the tree level correlator in the larger lattice size (L/a = 40, T /a = 80). In the long propagation length, the data can be fitted by a cosh function, f (x) = a cosh(b(x − 40)). At that time, the effective mass, which is shown in Fig. 2, is ω ∼ 0.144. In the short propagation length, the correlation function can be fitted by the power function, g(t) = c/t d . If we determine the fit parameter using the data at t/a = 2 and 4 by solving the equation, the exponent (d) becomes 3.59. There is a small discrepancy from d = 3 in Eq. (A4). We expect that it is an effect of the UV cutoff (=lattice spacing). Figure A2 is qualitatively consistent with the middle panel of Fig. A1. In our main analysis in this paper, we use the smaller lattice size than L/a = 40, T /a = 80, therefore the value of the effective mass (ω) is larger than the value of ω in L/a = 40. That is the reason why the shape of the correlation functions on our lattice shows cosh behavior, while we should stress that we take the continuum limit and the fermion mass is zero at the limit.

21/26
The other reason might come from the finiteness of the lattice extent. In two-dimensional (interactive) conformal field theory, the conformal map from the infinite plane to the cylinder with the radius (L) is known. On the cylinder coordinate, the correlation function becomes exponential function if the distant (in the direction with the infinite length of the cylinder) becomes larger than the compact radius L (See [34]), In our simulation, the temporal lattice extent is twice larger than the spatial one, and we found that the data in t/a ≥ L/2a regime shows the exponential behavior. Based on these considerations, we fit our lattice data using cosh function in Sec. 5 and Sec. 6.3.