RI/(S)MOM renormalizations of overlap quark bilinears with different levels of hypercubic smearing

On configurations with 2+1-flavor dynamical domain-wall fermions, we calculate the RI/(S)MOM renormalization constants (RC) of overlap quark bilinears. Hypercubic (HYP) smearing is used to construct the overlap Dirac operator. We investigate the possible effects of the smearing on discretization errors in the RCs by varying the level of smearing from 0 to 1 and 2. The lattice is of size $32^3\times64$ and with lattice spacing $1/a=2.383(9)$ GeV. The RCs in the $\overline{\rm MS}$ scheme at 2 GeV are given at the end, with the uncertainty of $Z_T$ reaching $\le1$% for the tensor current. Results of the renormalized quark masses and hadron matrix elements show that the renormalization procedure suppresses the $\sim$ 30% difference of the bare quantities with or without HYP smearing into the 3%-5% level.


Introduction
The uncertainties of renormalization constants (RCs) can be a main source of the uncertainties of hadronic matrix elements as lattice quantum chromodynamics (LQCD) is used to precisely calculate matrix elements numerically. For example, in our previous work on meson decay constants [1] the uncertainties of the tensor decay constants of D * (s) are dominated by the uncertainty of the RC of the tensor current (Z T ). Therefore it is crucial to compute RCs as precise as possible. In Ref. [2] the authors used momentumsubtraction schemes as intermediate schemes to determine Z T very precisely. After removing the nonperturbative effects from condensate contaminations they managed to obtain a precision better than 0.5% for Z T in the MS scheme at the scale of the bottom quark pole mass.
The momentum-subtraction scheme SMOM [3,4] was introduced to reduce the uncertainties of the RCs of quark bilinear operators, especially for the scalar density (Z S ) and therefore for the quark mass (Z m ) through Z m = 1/Z S for lattice chiral fermions. This scheme uses a symmetric combination for the momenta on the legs of the Green functions of the operators under renormalization. It is expected to have less nonperturbative infrared contaminations in the pseudoscalar and axial vector vertex functions compared with the MOM scheme [5] which uses forward Green functions. The perturbative conversion ratio of Z S from the SMOM scheme to the MS scheme was shown [6,7,8,9] to converge faster than that for the MOM scheme. Thus the truncation error from the perturbative conversion is reduced significantly.
In our previous work [10] on renormalization of quark bilinears for overlap fermions on domain-wall fermion configurations we compared SMOM and MOM schemes numerically. We indeed see that the infrared effects in the pseudoscalar and axial vector vertex functions are more suppressed in the SMOM scheme and Z S = Z P and Z V = Z A can be better satisfied. However we had difficulties to describe the a 2 p 2 dependence of Z MS S,P (2 GeV; a 2 p 2 ) obtained through the SMOM scheme, where a 2 p 2 is the initial scale in lattice units at which we do the renormalization and from which we run Z S,P to 2 GeV. It seems that we could not find a good renormalization window in a 2 p 2 in which both the nonperturbative effects and lattice discretization effects are small. While if we use the MOM scheme, a good window for Z MS S,P (2 GeV; a 2 p 2 ) can be found, in which the a 2 p 2 dependence can be well fitted by a linear function and is attributed to O(a 2 p 2 ) discretization effects. Similar observations were made in a more recent work [11], in which both the SMOM and MOM schemes were used to obtain the RCs for overlap quark bilinears on both domain-wall fermion and highly improved staggered quark configurations.
In the mixed action setup of the χQCD collaboration with overlap fermions on domain-wall fermion configurations, one level of hypercubic (HYP) smearing [12] on the gauge links is performed to construct the overlap Dirac operator. This can expedite the inversions of the Dirac operator [13]. Link smearing is widely used in LQCD simulations for various lattice fermions and has many benefits. For example, HYP smearing can suppress chiral symmetry breaking effects in Wilson fermions besides speeding up inversions for overlap fermions [14]. In Ref. [15] the authors studied the effects of smearing Table 1: Parameters of configurations used in this work. The residual mass of the domain-wall fermion in lattice units am res is in the two-flavor chiral limit as given in Ref. [18].
Label am l /am s Volume N conf am res f004 0.004/0.03 32 3  The upper end of the renormalization window is observed to be lowered by link smearing.
In this work we try to investigate the relation between the HYP smearing that we use to construct our overlap Dirac operator and the a 2 p 2 dependence of Z MS S (2 GeV; a 2 p 2 ) and also of other RCs. Both the MOM and SMOM schemes are used as intermediate schemes and the results are compared. The configurations used here were generated by the RBC/UKQCD collaborations with volume 32 3 ×64 and lattice spacing 1/a = 2.383 (9) GeV (labeled as 32I) [16]. This ensemble was also used in the study of [15]. We calculate the RCs for the tensor current and the quark field (Z q ) in the MS scheme on the 32I ensemble that were not done in our previous work [17]. The new three-loop conversion ratio for Z T from the SMOM scheme to the MS scheme [8] and four-loop ratios for Z S,T,q from MOM to MS are used in this work. Thus the systematic uncertainties of Z S,T are now reduced compared with those in Refs. [10,17].
The paper is organized as follows. In Sec. 2 we describe our lattice setup and briefly refresh the formulas of the MOM and SMOM schemes. Then in Sec. 3 we give the details of our calculation, the comparison of the results from different levels of HYP smearing and different intermediate schemes, and the discussions. We also investigate the renormalizated meson two-point functions at several quark masses in Sec. 4, to verify our renormalization procedure and estimate the discretization error due to different HYP smearing steps. Finally we summarize in Sec. 5

Simulation parameters
We use the 2+1-flavor gauge configurations generated by the RBC/UKQCD collaborations [16]. The sea quarks are domain-wall fermions with masses am l = 0.004, 0.006, 0.008 and am s = 0.03 for the mass degenerate up and down quarks and the strange quark, respectively. The parameters of the configurations are listed in Table 1. More detailed information can be found in Ref. [16].
For the valence quark we use overlap fermions. The massless overlap Dirac operator [19] is defined as where ε is the matrix sign function and D w (ρ) is the usual Wilson fermion operator, except with a negative mass parameter −ρ = 1/2κ − 4 in which κ c < κ < 0.25. We use κ = 0.2 in our calculation which corresponds to ρ = 1.5. The massive overlap Dirac operator is defined as To accommodate the SU(3) chiral transformation, it is usually convenient to use the chirally regulated fieldψ = (1 − 1 2 D ov )ψ in lieu of ψ in the interpolation field and the currents. This is equivalent to leaving the unmodified currents and instead adopting the effective propagator, where D c = ρDov 1−Dov/2 satisfies {γ 5 , D c } = 0 [20]. With the good chiral properties of overlap fermions, we can expect Z S = Z P and Z V = Z A . These relations were verified in our previous works on renormalizations of those bilinear quark operators [10,17].
We apply no smearing, one and two levels of HYP smearing in constructing the overlap Dirac operator as given in Eq.(1). In the rest of this paper the three cases are labeled as "thin," "HYP1," and "HYP2" respectively. Bootstrap and jackknife analyses are used to obtain the statistical uncertainties.
Same ten valence quark masses in lattice units are used for all smearing cases and ensembles. Their values are given in Table 2. The corresponding pion masses on the ensemble f004 are also given in Table 2, which are in the range of 220 to 500 MeV and are obtained by fitting two-point functions as described in Sec. 3.1. We use more statistics for the HYP1 case which corresponds to the setup used in most of the χQCD studies, and then the uncertainty of am π for this case is smaller. For the HYP1 case, we also list the m π L in the Table. The fitted pion masses at a same bare quark mass differ by renormalization effects and discretization errors.

Momentum subtraction schemes
The MOM scheme [5] is defined in the quark massless limit by imposing conditions on forward vertex functions Λ O of quark bilinear operators O at a renormalization scale µ where the subscript B stands for bare and Λ O can be calculated from the quark propagator S(p) and the Green function G O (p 1 , p 2 ) as with p 1 = p 2 = p. The projector Λ tree O (p) = Γ for the quark bilinear operatorsψΓψ (Γ = I, γ µ γ 5 , σ µν (= 1 2 [γ µ , γ ν ])) considered in this work. The trace "Tr" is over the color and spin indices. The quark field renormalization Z q is given by so that the chiral Ward identities are satisfied in the MOM scheme and we have Z MOM For the SMOM scheme [4] one considers symmetric momentum combinations in Eq. (5). The renormalization conditions for the scalar and tensor currents in the SMOM scheme are the same as Eq.(4) but at the symmetric point Eq. (7). The conditions for the quark field and the axial current are now and lim respectively. Here the subscript "sym" denotes the symmetric condition in Eq. (7). The conditions in Eqs. (8,9) are also consistent with the chiral Ward identities. Therefore one has Z SMOM A = Z MS A . The Green function G O is computed between two external off-shell quark states. Thus a gauge fixing (usually Landau gauge) is used to implement the renormalization conditions in both the MOM and SMOM schemes. In this study we use point source to calculate the quark propagator and then We use periodic boundary conditions in all four directions on our lattice when inverting the overlap Dirac operator. Thus the momentum modes in lattice units can be written as where L = 32, T = 64 and k µ are integers: For the MOM scheme we apply a democratic cut on the momentum modes to reduce the discretization effects from O(4) to H(4) symmetry breaking on the lattice. Note this cut is more strict than those used in Refs. [10,17]. It is difficult to apply the above democratic cut on all the symmetric momentum modes (q 2 = p 2 1 = p 2 2 ) required by the SMOM scheme. But still we can apply a cut on both p 1 and p 2 p [4] (p 2 ) 2 < 0.34, for p 1 and p 2 .
That is to say p 1 and p 2 are almost body diagonal and q is almost along an axis (to satisfy q 2 = p 2 1 ). A typical example (with p [4] /(p 2 ) 2 = 0.25) is The cut 0.34 in Eq. (14) is looser than the 0.26 in Eq.(13) because the symmetric condition q 2 = p 2 1 should also be satisfied. If 0.26 is applied, too few momentum modes will be left. In this way the numerical results of RCs from the SMOM scheme will be less scattered around a smooth curve.

Local axial vector current
We start with the renormalization of the local axial vector current. The RC Z A is then used to scale the other RCs for the quark field, the scalar density and the tensor current. Z A can be obtained by using the partially conserved axial current relation where A µ and P are the axial vector current and pseudoscalar density, respectively, and the quark mass RC Z m equals Z −1 P for lattice chiral fermions such as the overlap fermion used in this work.
Consider the zero momentum two-point correlators in the pseudoscalar channel From Eq.(16) one finds and or where the partial derivative is replaced by a difference. One can simultaneously fit the two-point correlation functions in Eqs. (17,18) at large time, where the ground state contribution dominates, and use Eq. (19) to calculate Z WI A . Alternatively one may use the ratio Eq.(21) to get Z WI A . We compared the results from the two methods in our previous work [10] and found that Z A from the first method has less dependence on the valence quark mass. Thus, we choose to use the first method in this work.
Two-point correlation functions on 42 configurations are calculated with a Z3 noise wall source at a given time slice for ensemble f004 for the three smearing cases. Note no gauge fixing is needed for calculating the gauge invariant correlation functions.
The results of Z WI A are collected in table 3 and drawn in figure 1. To fairly compare the three smearing cases we use the same statistics in figure 1. We can clearly see that smearing decreases the statistical error of Z WI A and drives it closer to one as the level of smearing increases. By linearly extrapolating Z WI A to the valence quark massless limit, one obtains Z WI A = 1.445(13), 1.0808(10), 1.0506(6), for thin, HYP1 and HYP2, respectively. Since the uncertainty of Z WI A using the thin link is much larger than the other cases, we repeat the calculation with a Coulomb wall source [11] to improve the statistics, and obtain Z WI A,thin = 1.4403(6) which is consistent with that listed in Eq. (22) but has much smaller uncertainty.
We also compare Z WI A for the case HYP1 obtained here with that from our previous work in figure 1. One can see that the results are in good consistency although different boundary conditions in the time direction are used here and in [17]. This is not surprising since RCs are not sensitive to finite volume effects from different boundary conditions.  (9) 1.0804 (9) 1.0803(6) HYP2 1.0506 (7) 1.0507 (7) 1.0507 (6) 1.0507 (5) 1.0509(4)  ) for all three smearing cases. The straight lines in both graphs are linear chiral extrapolations of Z A . The values in the chiral limit are given in the right graph. For the case HYP1 comparison is also shown with our previous work [17] with an antiperiodic boundary condition in the time direction.
We will need the RCs for the HYP1 case for calculating hadronic matrix elements in the future. To shrink the statistical errors for this smearing case we use 628 configurations, 16 sources on 42 and 16 sources on 49 configurations to compute the two-point correlators in Eqs. (17,18) on ensembles f004, f006 and f008, respectively. The linear chiral extrapolation of Z WI A (am q ) is similar to those shown in figure 1. In the chiral limit of valence quark we find respectively, on the three ensembles. Here the error is statistical and we used all ten valence quark masses for doing the chiral extrapolations. For the chiral limit of the up and down quarks in the sea we linearly extrapolate the results on the three ensembles to am l + am res = 0 and obtain Z A = 1.0789 (7). This extrapolation as well as a constant fit (1.0788 (2)) is shown in figure 2. The difference between the linear extrapolation and the constant fit is taken as a systematic error (0.0001) below. We repeat the above chiral Another systematic uncertainty comes from the massive strange quark in the sea. By using half of the slope from the linear chiral extrapolation of the up/down sea quarks, we estimate the change of Z A to be 0.0003 in the limit am s + am res = 0. In the end we have Z A = 1.0789(7)(1)(6)(3). Adding up all the errors quadratically, one gets Z A = 1.0789 (10). This value agrees with our previous result 1.086(2) [17] at ∼ 3σ. We now have more statistics and consider more carefully the systematical uncertainties in our new result.
From the fitting of the two-point functions we also obtain the pion mass at each valence quark mass as given in table 2.
Since the renormalization conditions in both the MOM and SMOM schemes are consistent with the chiral Ward identities, we have Z MOM A in the continuum limit. On the lattice they can differ by discretization effects.

Scalar density
The scalar density renormalization constant Z MOM S can be obtained from Z WI A computed in Sec. 3.1 and the ratio of projected vertex functions, where The scale dependence of the ratio in Eq. (24) is governed by the anomalous dimension of the scalar density since Z A is scale independent. The chiral extrapolation of the valence quark is done with an ansatz which was also used in Refs. [3,10,11,17]. B s is taken as the chiral limit value of Z S /Z A . We find the contribution of the first term on the right-hand side of Eq. (26) in our data is small. Thus, we also tried linear extrapolations in am q and checked that consistent results are obtained in the chiral limit. Examples of both extrapolations at a 2 p 2 = 3.855 or 8.135 on ensemble f004 are shown in figure 3. In our following analyses we generally use data points at scales above a 2 p 2 = 4 to avoid possible large nonperturbative effects. For the three smearing cases we obtain Z MOM S /Z MOM A in the valence quark chiral limit by using Eq. (26). The results are shown in the left panel of figure 4.
The conversion to the MS scheme is done by using the ratio in Landau gauge up to four loops [21,22,23] where n f is the number of flavors. To run the results in the MS scheme from the initial scale |p| to 2 GeV, we use the quark mass anomalous dimension given in Ref. [22] since The four-loop (next-to-next-to-next-to-leading order) running results are shown in the right panel of figure 4. Here we use the inverse lattice spacing 1/a = 2.383 (9) GeV as determined in [24] to get the position of 2 GeV. This 1/a is about 4% higher than that used in our previous work [17]. Λ MS QCD = 332 (17) MeV [25] is used to calculate the strong coupling constant by using its perturbative running to four loops [26,27]. Now we turn to the SMOM scheme. The ratio for the three smearing cases is shown in the left panel of figure 5 in the valence quark massless limit. Here The massless limit is obtained from linear extrapolations in the valence quark mass. The perturbative conversion ratio to the MS scheme has been calculated up to three loops [6,7,8,9] (88) After the conversion we can run Z MS S (a 2 p 2 )/Z A to 2 GeV. The results are shown in the right panel of figure 5.

Comparison of the three smearing cases (via MOM)
As shown in the right panel of figure 4, the dependence of Z MS S (2 GeV, a 2 p 2 )/Z A on a 2 p 2 is similar for all three cases when MOM is used as the intermediate scheme. At large momentum scale this a 2 p 2 dependence can be described by a polynomial function where the last two terms on the right-hand side contain the discretization effects in Z MS S (2 GeV, a 2 p 2 )/Z A . The fittings of Eq. (31) without the c 2 term to data points in the range a 2 p 2 ∈ [5, 10] give χ 2 /dof < 1 for all three smearing cases. Therefore the O((a 2 p 2 ) 2 ) discretization effects are smaller than our current statistical error. We take this linear extrapolation value Z MS S (2 GeV)/Z A as the final result, which is collected in table 4. The fitting range of a 2 p 2 is varied to [4,10] and [6,10] to estimate the systematic error in Z MS S (2 GeV)/Z A given in table 6 from the choice of the lower bound. There are curvatures in Z MS S (2 GeV, a 2 p 2 )/Z A for HYP1 and HYP2 especially at a 2 p 2 < 5. Thus we also try the following ansatz to study the a 2 p 2 dependence for these two smearing cases We find both Eq. (31) and Eq. (32) can describe the data well in the range a 2 p 2 ∈ [2, 10] with a roughly same χ 2 /dof, which is less than 1. Thus, it is hard to attribute the curvature to only nonperturbative effect or O(a 2 p 2 ) effect. We then use Eq. (33) to fit the data in a 2 p 2 ∈ [5, 10]. The resulted center values of Z MS S (2 GeV)/Z A change by 0.82% and 0.41% for HYP1 and HYP2, respectively. This ansatz dependence is checked on all three ensembles f004, f006 and f008 for HYP1. The resulted change (0.98%) in the final result of Z MS S (2 GeV)/Z A in the sea quark massless limit is given in table 6 as a systematic error.

Comparison of the three smearing cases (via SMOM)
When the SMOM scheme is used as the intermediate scheme, the dependence of Z MS S (2 GeV, a 2 p 2 )/Z A on a 2 p 2 (right panel of figure 5) looks similar for the two smearing cases HYP1 and HYP2, which, however, seems to be quite different from the thin link case. The a 2 p 2 dependence is more flat at a 2 p 2 > 4 when no smearing is applied. All linear extrapolations in a 2 p 2 in the ranges a 2 p 2 ∈ [4, 10], [5,10] and [6,10] have χ 2 /dof < 1 and give consistent extrapolated results for the thin link case. The numbers from using the range [5,10] are given in table 4.
If we use Eq. (31) to fit Z MS S (2 GeV, a 2 p 2 )/Z A with the c 2 term for the HYP1 case, the χ 2 /dof are, respectively, 1.78, 0.78 and 1.06 for a 2 p 2 ∈ [4, 10], [5,10] and [6,10]. If we drop the c 2 term, that is to say, use a linear fit, then we get large χ 2 /dof (42, 8.7 and 3.7 for the three ranges, respectively). Thus we choose the second order polynomial extrapolation and get Z MS S (2 GeV)/Z A = 0.9126 (37) in table 4. The dependence on extrapolation ansatz is checked by using Eqs. (32,33) on all three ensembles f004, f006 and f008. The variation in the center value is treated as a systematic error as collected in table 6.
For the case HYP2 of SMOM we tried six ansatz: linear, second and third order polynomials of a 2 p 2 with or without an inverse term c −1 /(a 2 p 2 ) to take into account possible nonperturbative effects at small a 2 p 2 . None of them can easily fit the data in the above three ranges of a 2 p 2 . The χ 2 /dof are in the range 2.3 to 107.
One possible reason is the following. The statistical error decreases as the level of smearing increases. Then the other systematic errors (such as the O(4) to H(4) breaking effects) are no longer small compared with the statistical error. We use the (near) diagonal momentum modes to suppress those effects. However, in the case of SMOM more off-diagonal momentum modes have to be used than in the MOM case. Therefore, we increase the error of the fitting results by χ 2 /dof to compensate the ignorance of those systematic effects.
If we use the same range a 2 p 2 ∈ [3. 5,9] as that in Ref. [11] for ensemble 64I and use a same third order polynomial, then we obtain χ 2 /dof = 41.4/17 = 2.4 and Z MS S (2 GeV)/Z A = 0.9657(88), where the error has been enlarged by χ 2 /dof. The a 2 p 2 dependence of Z MS S (2 GeV, a 2 p 2 )/Z A obtained through the SMOM scheme for the case HYP2 cannot be readily described by polynomial discretization effects even at large a 2 p 2 . This dependence seems to change significantly as one applies smearing compared with the thin link case. For the case HYP1 the extrapolated results are more sensitive to the fitting range of a 2 p 2 as shown by the corresponding systematic uncertainty in table 6 compared with MOM. When MOM is used as the intermediate scheme, the a 2 p 2 dependence of Z MS S (2 GeV, a 2 p 2 )/Z A can be well described by a linear function for all three smearing cases at large a 2 p 2 . Thus we use the results obtained Table 5: Z MS S (2 GeV)/Z A on ensembles f004, f006 and f008 for HYP1. They are obtained after removing the discretization effects by using a straight-line (via MOM) or a secondorder polynomial (via SMOM) extrapolation in a 2 p 2 . The error is from statistics and a 2 p 2 extrapolation. The values in the sea quark massless limit are from linear extrapolations. For the intermediate MOM scheme the error is enlarged by

Final results for one level of HYP smearing
On the other two ensembles f006 and f008 with smearing HYP1, we do similar calculations in the MOM scheme and obtain Z MS S (2 GeV)/Z A = 0.9287 (26) and 0.9352 (17), respectively, in the valence quark massless limit. The results on all three ensembles are collected in table 5 as well as those obtained through the SMOM scheme. Then with the results on all three ensembles we do a linear extrapolation in (am l + am res ) to the light sea quark massless limit and find Z MS S (2 GeV)/Z A = 0.9348(88) with a confidence level 0.021, where the error has been enlarged by χ 2 /dof = √ 5.3. The slope from this extrapolation is consistent with zero (−0.09(57)), of which the center value is used to estimate the systematic uncertainty due to a nonzero strange sea quark mass given in table 6. We also tried a constant fit of Z MS S (2 GeV)/Z A on the three ensembles to go to the sea quark massless limit. The difference in the center values from the linear and constant fits is set to be the systematic error for extrapolation in m l in table 6.
We list all systematic uncertainties in table 6 for Z MS S (2 GeV)/Z A . They are estimated in similar ways as in Refs. [17,10,11]. The truncation error in the conversion ratio from MOM to MS is found to be 0.67%, which is smaller than the 1.5% in [17] and the 2.29% in [11]. One reason is that here the lower limit of |p| used in the a 2 p 2 extrapolation is about 5.3 GeV, which is higher than the 4 GeV in [17] and the 3 GeV in [11]. The other reason is now we use the newly calculated four-loop conversion ratio Eq. (27). Each of Λ MS QCD and the inverse lattice spacing is varied in one sigma to check the resulted change in Z MS S (2 GeV)/Z A . The perturbative running in the MS scheme includes a four-loop result and thus gives a negligible truncation error. In total the systematic uncertainty is found to be 1.28% when one uses the MOM scheme.
By using Z A = 1.0789(10) from section 3.1, we then find Z MS S (2 GeV) = 1.009(16). This center value is about 4% (or 1.6σ) away from our previous result 1.056 (6)(24) given in [17]. This change is mainly from the different inverse lattice spacings used here and there. Also, the a 2 p 2 extrapolation ranges are different in the two works. We confirmed that we obtain consistent results if the same inverse lattice spacing and a 2 p 2 range are used here as in [17]. Through the SMOM scheme one gets Z MS S (2 GeV) = 0.988(53), which agrees with that obtained through the MOM scheme. Our final result Z MS S (2 GeV) = 1.009(16) agrees with the 1.034 (25) given in [11] for its ensemble 64I, which has a similar setup as the ensemble used here.

Tensor current
The renormalization constant of the tensor current is needed in calculating observables such as the tensor decay constant of vector mesons [1] and nucleon isovector tensor charge [28]. The ratio of the RCs of the tensor and axial vector current in the MOM scheme is given by where Γ T (p) = 1 144 Tr The left panel of figure 6 shows the numerical results of the above ratio for the three smearing cases on ensemble f004 in the linearly extrapolated valence quark massless limit.
The conversion ratio to the MS scheme for Z MOM T up to four loops is given by [29, After finishing the conversion and four-loop running to 2 GeV from the initial scale a 2 p 2 in the MS scheme, we obtain the right panel of figure 6. The calculation in the SMOM scheme starts with the ratio where The numerical results of this ratio are shown in the left panel of figure 7 for all three smearing cases in the valence quark massless limit on ensemble f004. The valence chiral extrapolation is done by using a linear function in am q .
By using this conversion ratio and the anomalous dimension of Z T up to four loops in the MS scheme [31] we obtain Z MS T (2 GeV; a 2 p 2 )/Z A as a function of the initial scale a 2 p 2 . The results are shown in the right panel of figure 7.
For n f = 3 and scale |p| = 5.3 GeV, the conversion ratio in Eq. (40) is Assuming the coefficient in the O(α 4 s ) term is 3 times as big as that in the O(α 3 s ) term, we can estimate the size of the O(α 4 s ) term to be around 0.0026. Therefore, the truncation error is about 0.26%, which is a little bit larger than that for the MOM scheme at the same order. One can then extrapolate Z MS T (2 GeV; a 2 p 2 )/Z A at large a 2 p 2 to a 2 p 2 = 0 to remove lattice artifacts proportional to a 2 p 2 (and (a 2 p 2 ) 2 ) by using a linear function (or a secondorder polynomial) of a 2 p 2 . To avoid the nonperturbative effects at small a 2 p 2 we use Z MS T (2 GeV; a 2 p 2 )/Z A in the range a 2 p 2 ∈ [5,10]. We vary the range to [4,10] and [6,10] to estimate the systematic error from the choice of fitting ranges.
For the case HYP1 we find that a linear function of a 2 p 2 is good enough (χ 2 /dof < 1) to describe Z MS T (2 GeV; a 2 p 2 )/Z A in the range a 2 p 2 ∈ [5,10] or [6,10] obtained through either the MOM or the SMOM scheme. Figure 8 shows the linear extrapolations using the data at a 2 p 2 ∈  in the range a 2 p 2 ∈ [4, 10] we get χ 2 /dof = 1.7, which is more or less acceptable, and 0.7 for MOM and SMOM, respectively. The extrapolated results are given in table 7 along with the χ 2 /dof for different fitting ranges. The uncertainties in the table are from statistics and the extrapolations. The change in Z MS T (2 GeV)/Z A is around or less than 0.2% as we vary the fitting ranges for both intermediate schemes. The difference in Z MS T (2 GeV)/Z A from the two schemes is around 0.9%. For the case without smearing, the linear extrapolations of Z MS T (2 GeV; a 2 p 2 )/Z A to a 2 p 2 = 0 are similar to those for HYP1. The χ 2 /dof are smaller because the statistical uncertainties in the data are bigger. The difference in Z MS T (2 GeV)/Z A from the two schemes is around 1.7%.
For the case HYP2 we find that linear extrapolations in a 2 p 2 give large χ 2 /dof (> 2) when MOM is used as the intermediate scheme. Thus a second-order polynomial in a 2 p 2 is used for the extrapolation (χ 2 /dof = 16/32 for a 2 p 2 ∈ [5, 10]) and we obtain  Table 8: Z MS T (2 GeV)/Z A on ensembles f004, f006 and f008. They are obtained after removing the discretization effects by using a straight-line extrapolation in a 2 p 2 . The error is from statistics and a 2 p 2 extrapolation. The values in the sea quark massless limit are from linear extrapolations, whose errors are enlarged by Z MS T (2 GeV)/Z A = 1.0559 (7). For the intermediate SMOM scheme both linear and second-order polynomial extrapolations have χ 2 /dof ∼ 1.6 and give consistent extrapolated results 1.0705(9) and 1.0731(54). The difference in Z MS T (2 GeV)/Z A from the two schemes is 1.6% from the second-order polynomial extrapolation.
We do similar analyses of Z MS T (2 GeV)/Z A for the other two ensembles f006 and f008 for the case HYP1. The final results for all three ensembles are collected in table 8. The values in the sea quark massless limit are from linear extrapolations in am l + am res . The two numbers from the two intermediate schemes are in agreement at 1.2σ, with a 1.0% change in the center values.
The systematic uncertainties of Z MS T (2 GeV)/Z A are given in table 9. The analysis procedure is similar to that for Z MS S (2 GeV)/Z A . To check the model dependence of the a 2 p 2 extrapolation we also use the following ansatz besides the linear extrapolation. The resulted difference in the chiral limit is included where the two errors are statistical and systematic, respectively. We take the number which is in agreement with the result 1.150(5) for ensemble 64I in Ref. [11] (table XVII).

Quark field renormalization
QCD propagators such as the quark propagator can provide nonperturbative information about QCD as, for example, discussed in Ref. [32]. Studies of quark propagators in Landau gauge by using lattice QCD can be found in, for example, Refs. [33,34]. The quark field renormalization constant Z MS q in the MS scheme is useful when quark propagators are used to determine the quark chiral condensate as was tried in Ref. [35].
To obtain Z MS q the following ratios are calculated and then converted to the MS scheme: Figure 9: Extrapolations of Z MS q (2 GeV; a 2 p 2 )/Z A to a 2 p 2 = 0 from the two intermediate schemes. a 2 p 2 is the initial renormalization scale squared.
From the linear or second-order polynomial dependence on a 2 p 2 at large scale (a 2 p 2 ∈ [5, 10]) we extrapolate Z MS q /Z A (2 GeV; a 2 p 2 ) to a 2 p 2 = 0 to remove the lattice artefacts. Besides the linear function or Eq. (50), we try functions with a nonperturbative term to estimate the uncertainty of Z MS q (2 GeV)/Z A from the fitting ansatz as given in table 11.
The extrapolated results Z MS q /Z A (2 GeV) on the three ensembles with HYP1 are listed in table 10 along with the values in the sea quark massless limit. The sea quark chiral extrapolation is done linearly in (am l + am res ). A constant extrapolation in the light sea quark mass is used to estimate the associated systematic error. The slope from the linear extrapolation is used to estimate the systematic error due to the nonzero strange sea quark mass. The final results in the chiral limit obtained through the two intermediate schemes are in good agreement.
The systematic uncertainties of Z MS q (2 GeV)/Z A are given in table 11. Then we find where the first error is statistical and the second error systematic. The dominant error is statistical due to the a 2 p 2 -extrapolation with a second-order polynomial. Adding up all the errors quadratically and using Z A = 1.0789(10) from section 3.1, we get Table 10: Z MS q (2 GeV)/Z A on ensembles f004, f006 and f008 for HYP1. They are obtained after removing the discretization effects by using a second-order polynomial function in a 2 p 2 . The error is from statistics and a 2 p 2 extrapolation. The values in the sea quark massless limit are from linear extrapolations, whose error is enlarged by

Test the renormalization
In this section, we check the HYP smearing dependence of the renormalized light quark mass and hadron matrix elements on the f004 ensemble, based on the renormalization constants obtained in the previous section and collected in Table 12 (without the chiral extrapolation of the sea quark mass). Besides the systematic uncertainty from the sea quark mass extrapolation, we also dropped the systematic uncertainty from the conversion ratio for Z S since it is independent of the HYP smearing setup. The quark mass dependence of the pion mass squared is shown in figure 10 on the ensemble f004. About 42 configurations are used for all three smearing cases. In the left panel we do linear fits (am π ) 2 = A · am q + B separately for the three cases, where m q is the bare quark mass. The data can be well described by this linear function and the intercepts B are all consistent with zero with our statistical uncertainties. At a same bare quark mass the pion masses are different because the renormalized quark masses are different (and the discretizaton effects are also different) for the different smearing cases.
In the right panel of figure 10 we plot the pion mass squared as a function of the renormalized quark mass m R q (2 GeV; MS) = m q /Z MS S (2 GeV), where Z MS S (2 GeV) are given in table 12. Then we find that all data points can be fitted to one linear function (am π ) 2 = A ′ · am R q + B ′ . Thus, after applying renormalization for the different smearing cases one can expect to obtain a renormalized light quark mass independent of smearing with the given statistics.
Because of the correlation between the data with different HYP smearing steps, the independence of smearing can be checked with even higher precision. We calculate the point 4-4-4 grid source [13] propagators on two time slices on about 40 configurations, and construct the meson correlators with different gamma matrices. Such a setup is equivalent to 4 3 × 2 = 128 measurements per configuration of the following point source where m i and ⟨O|i⟩ are the ith state mass and matrix elements, respectively. Since the HYP smearing can change the UV behavior of the hadron spectrum, we will concentrate on the impact of HYP smearings on the ground state mass m 0 and also the renormalized ground state matrix element ⟨O|0⟩ R = Z O ⟨O|0⟩. On a lattice with finite size T in the time direction, Eq. (55) should be modified into For the pseudoscalar correlator which has good signal around T /2, we apply the one-state fit with the following ansatz in the range 0 ≪ t ≪ T , and do the two-state fit for the cases of other hadrons, at relatively smaller t as the results around T /2 can be very noisy, where c 1 and δm are additional parameters to describe the contaminations from the excited states. We do the foldingC(t) = 1 2 [C(t) + C(T − t)] on the correlator and require theC(t) at the maximum t used in the fit to have at least 3-5σ signal (or T /2 in the pseudoscalar case), and tune the minimum t to make the χ 2 /d.o.f. to be around one and the corresponding Q value of the fit to be larger than 0.05.
Because of the discretization error, the meson mass using the same renormalized quark mass can be different with different HYP-smearing steps. Thus we can consider this effect in a reversed way, by tuning the quark masses to make the corresponding pseudoscalar masses obtained by C 2,P as the following four values regardless of the HYP-smearing steps: 1) 0.302 GeV which corresponds to the unitary pion mass on the f004 ensemble; 2) 0.675 GeV which corresponds to the strange quark mass, m MS q (2 GeV) ∼ 0.1 GeV; 3) 0.976 GeV which corresponds to m MS q (2 GeV) ∼ 0.2 GeV; 4) 1.230 GeV which corresponds to m MS q (2 GeV) ∼ 0.3 GeV. Practically we calculate at two quark masses around each of the above cases, and do the interpolation to make the pseudoscalar mass to be exact.
Then we compare the renormalized quark mass with different HYP-smearing steps. The ratio of the quark masses with one or two steps HYP smearing over that with the thin link are plotted in Fig. 11. The results with different HYP smearing are based on the same configurations, and the data correlation has been taken into account to suppress the uncertainty of the ratios. of the renormalized quark masses with n = 1 (yellow triangles) and n = 2 (green dots) steps of HYP smearing, to make the corresponding pseudoscalar meson mass to be 0.302, 0.675, 0.976 and 1.23 GeV, respectively.
In Fig. 11, the yellow triangles show the ratio of the quark mass with one-step HYP smearing (the standard χQCD setup). We can see that the ratio is around 2%-3% smaller than 1. Such a difference is slightly larger than the joint uncertainty of Z S with the thin link, but an order of magnitude smaller than that of the bare quark masses with or without HYP smearing (∼ 20%). On the other hand, the green dots show the case using the quark mass with two-step HYP smearing in the numerator. The deviation seems to be smaller comparing to the one-HYP case, and the result at unitary pion mass seems to be consistent with one with much larger statistical uncertainty. Thus, we conclude that the HYP smearing may introduce around 3% discretization error in the renormalized quark mass at this specific lattice spacing.
Then we are ready to compare the renormalized ground state matrix elements ⟨O|0⟩ with different O, through the ratio R n=1,2 O = ⟨O|0⟩ nHYP R /⟨O|0⟩ thin R . R = 1 means the renormalization matrix elements are independent of the HYP smearing steps, and a nonvanishing deviation suggests a discretization error as the HYP smearing dependence should vanish in the continuum.
for the pseudoscalar (left panel, O =qγ 5 q) and scalar (right panel, O =qq) matrix elements using Z S = Z P . The lightest quark mass case in the scalar channel is dropped due to the mixture between the single-hadron and multiple-hadron states. Fig. 12 shows the cases of the pseudoscalar (left panel) and scalar (right panel) matrix elements, while the pseudoscalar case has much better signal and shows an obvious quark mass dependence. After the chiral extrapolation, the deviation from one in the one-HYP pseudoscalar case is around 3%, while the two-HYP case is around 1%. The uncertainty in the scalar case is much larger and the result after chiral extrapolation is consistent with the pseudoscalar case. Note that the lightest quark mass case in the scalar channel is dropped since the mixture between the σ and ππ states makes the two-state fit defined in Eq. (58) to be unreliable. For similar reason, we dropped the light quark mass cases in all the channels except the pseudoscalar ones.
In the channels using the vector or axial-vector currents, the pseudoscalar channel using O =qγ 5 γ 4 q also has good signal (upper left panel of Fig. 13) and suggests that the renormalized matrix elements in the chiral limit is independent of the HYP smearing, for both the one-HYP and two-HYP cases. Similarly to the scalar case, the results from the spatial components of the vector and axial-vector currents are noisier as shown in the other panels of Fig. 13, while the discretization error due to the HYP smearing is at the 5% level in the cases we investigated.
As shown in Fig. 14, the operator O =qγ 4 γ i q corresponds to the vector channel which is lighter than the tensor channel using the operator O =qγ i γ j q, and then has a better signal. We can see that the HYP smearing dependence in the vector channel is at the same level as the other cases, while the tensor channel seems to have larger deviations but the uncertainties are also large.

Summary
Renormalization constants are necessary in lattice QCD calculations of various hadronic matrix elements, which are important in precise determinations of the parameters of the Standard Model and in searching new physics. Thus, one needs to calculate RCs as precise as possible.
Gauge link smearings are widely used in lattice QCD calculations. They can suppress the ultraviolet fluctuations of the gauge fields and decrease the statistical uncertainties in practice calculations besides bringing in many other benefits. They can also affect the vertex functions of quark operators which are used in (S)MOM renormalization of those operators. Therefore we investigate the effects of HYP smearing on the renormalization window in this work as we try to get RCs as precise as possible by using both MOM and SMOM as intermediate schemes. We check the a 2 p 2 dependence of Z MS X (2 GeV; a 2 p 2 ) (X = S, T, q) as it is obtained by running from the initial scale |p| to 2 GeV.
For our lattice setup, in general, the a 2 p 2 dependence of Z MS X (2 GeV; a 2 p 2 ) via the MOM scheme is not very sensitive to the HYP smearing although it is affected by the levels of smearing that we use (no smearing, one hit and two hits). The a 2 p 2 dependence of Z MS S,T (2 GeV; a 2 p 2 ) can be described straightforwardly by a linear function in the range a 2 p 2 ∈ [5, 10] for all smearing cases. In physical units the range is |p| ∈ [5.3, 7.5] GeV. The a 2 p 2 dependence of Z MS q (2 GeV; a 2 p 2 ) can be described by a second-order polynomial for all three smearing cases.
For the intermediate SMOM scheme, the behaviour of the a 2 p 2 dependence of Z MS X (2 GeV; a 2 p 2 ) is sensitive to whether HYP smearing is applied or not. From the right panels of figures 5, 7 and 9 we can see that the behaviors of Z MS X (2 GeV; a 2 p 2 ) for HYP1 and HYP2 are apparently different from that for the thin case. Z MS X (2 GeV; a 2 p 2 ) can be described by a linear function of a 2 p 2 when no smearing is applied. Among the three RCs (X = S, T, q) the a 2 p 2 dependence of Z MS T (2 GeV; a 2 p 2 ) has the least sensitivity to smearing. It can still be fitted to linear functions for both HYP1 and HYP2 although the slope has a sizable change compared with no smearing. The a 2 p 2 dependence of Z MS q (2 GeV; a 2 p 2 ) changes from a linear function to a second-order polynomial after smearing is applied. Z MS S (2 GeV; a 2 p 2 ) has the largest sensitivity to smearing. With HYP2 we find it is hard to describe the a 2 p 2 dependence of Z MS S (2 GeV; a 2 p 2 ) no matter if we use high order terms for discretization effects or inverse terms for possible nonperturbative effects.
By going to higher renormalization scale than in [17] and using four-loop conversion ratios from the MOM scheme to MS, we reduce the systematic errors of the RCs for the scalar and tensor operators. For Z MS T (2 GeV) we obtain a total uncertainty less than 1%.
We also checked the HYP smearing dependence of the renormalized quark masses and hadron matrix elements. The results show that the renormalization suppresses the ∼ 30% difference of the bare quantities into 3%-5% level, while more than one step of HYP smearing may not make the residual deviation to be smaller.