$O(4)$-symmetric position-space renormalization of lattice operators

We extend the position-space renormalization procedure, where renormalization factors are calculated from Green's functions in position space, by introducing a technique to take the average of Green's functions over spheres. In addition to reducing discretization errors, this technique enables the resulting position-space correlators to be evaluated at any physical distance, making them continuous functions similar to the $O(4)$-symmetric position-space Green's functions in the continuum theory but with a residual dependence on a regularization parameter, the lattice spacing $a$. We can then take the continuum limit of these renormalized quantities calculated at the same physical renormalization scale $|x|$ and investigate the resulting $|x|$-dependence to identify the appropriate renormalization window. As a numerical test of the spherical averaging technique, we determine the renormalized light and strange quark masses by renormalizing the scalar current. We see a substantial reduction of discretization effects on the scalar current correlator and an enhancement of the renormalization window. The numerical simulation is carried out with $2+1$-flavor domain-wall fermions at three lattice cutoffs in the range 1.79--3.15~GeV.


I. INTRODUCTION
Operator renormalization is necessary to calculate many quantities such as weak matrix elements using lattice simulation. So far, several different methods to renormalize lattice operators have been proposed, applied and improved. Since each method has individual advantages and disadvantages, we can use the method that is the most convenient for our situation and purpose.
In this work, we focus on the position-space procedure [1,2], in which the renormalization condition for an operator is imposed on the corresponding correlation function in position space. It is an important advantage of this procedure that it provides a fully gauge invariant renormalization prescription since the correlator used in the renormalization condition is gauge invariant. This advantage prevents the mixing with gauge noninvariant operators that occurs in gauge noninvariant schemes such as the regularization independent momentum subtraction (RI/MOM) scheme [3]. Since the operators appearing in the position-space renormalization prescription are evaluated at separated space-time points, operators which vanish when the equations of motion are imposed will also not contribute. Therefore a position-space renormalization scheme will also avoid mixing with operators which vanish by the equations of motion -mixing which can occur in the RI/MOM approach.
An important difficulty of the position-space approach arises from the discrete lattice of points on which the position space Green's function is evaluated. Unless one works with lattices whose lattices spacings are related as integer multiples, errors may be introduced when combining results from two different ensembles. Combining results from ensembles with different lattice spacing is necessary both when evaluating the continuum limit and when using step scaling [4][5][6]. (For example, when Cichy et al. [7] employ step scaling in position space they consider lattice spacings which differ by factors of two.) Recall that step scaling is an important nonperturbative method used to relate the normalization of operators that are being used in a coarse lattice calculation to physically equivalent operators defined on a fine, weak-coupling lattice where a connection to perturbatively normalized operators can be more accurately made.
In an RI/MOM scheme the Fourier transform averages over the discrete lattice and the resulting functions of momentum approach their continuum limits in a well-understood fashion [8,9]. In this paper we propose an alternative average that partially smooths the discrete nature of the position-space lattice while working with gauge invariant quantities and maintaining a non-zero separation between the operators whose Green's functions are being studied.
Our strategy is best illustrated using two-point functions, which are the starting point of the present position-space renormalization schemes, as illustrated by the expression where O is a gauge-invariant composite local operator, x n is a point on our discrete lattice determined by the four integers n = (n 1 , n 2 , n 3 , n 4 ) and, for simplicity, the second point 0 is chosen to be the origin, also a point of this lattice. As is described in more detail in Section III, we begin by extending this function into a function G(x) of the continuous position four-vector x, obtained by multi-linear interpolation from the sixteen values obtained by evaluating G(x n ) at the sixteen lattice points that lie at the vertices of the four-dimensional cube which contains the point x.
Assuming, as we do throughout this paper, that our lattice theory has no order a errors, our interpolated Green's function G(x) will agree with the corresponding Green's function of the continuum theory up to errors which vanish as a 2 in the continuum limit. Of course, the a 2 errors which appear in the piece-wise linear function G(x) will still reflect the O(4) symmetry breaking of the underlying lattice. In order to reduce these lattice artifacts and define a function of a single scale, we further simplify our Green's function by averaging the point x over a three-dimension sphere of radius |x| centered at the origin: We will then impose conditions on G(|x|) to renormalize the operator O.
The finite lattice spacing errors that are present in the lattice Green's function G(x n ) are expected to appear as simple polynomials in a with coefficients which in perturbation theory depend only logarithmically in a, allowing a simple exptrapolation to the continuum limit.
The lattice spacing dependence of our averaged quantity G(|x|) will be more complicated.
In addition to the simple a 2 errors coming from G(x n ), the sphere averaging procedure will introduce O(a 2 ) errors which, while bounded by a 2 may be complicated irregular functions of a which could cause an explicit extrapolation in a 2 to fail. As will be shown in Appendix A, this irregular dependence on a appears to be negligible, making the scheme proposed here suitable for a calculation in which the continuum limit is to be evaluated. Of course, were this error too large, we may be able to introduce a higher-order interpolation scheme which would make these troubling effects of higher order than a 2 and therefore systematically negligible.
As an example of the spherical average, we present our result for the quark mass renormalization, which can be done by renormalizing the scalar current. There are several previous works on position-space renormalization of bilinear operators [2,10,11]. For renormalization of bilinear operators, there is another important advantage of the position-space procedure: the perturbative matching to the modified minimal subtraction (MS) scheme is available to O(α 4 s ) for the vector, axial-vector, scalar and pseudoscalar currents and to O(α 3 s ) for the tensor current [12]. Utilizing the spherical averaging technique, we perform a new analysis that takes the continuum limit of the renormalized quark mass at many values of |x| and shows its |x|-dependence. The final result agrees with the FLAG average [13] as well as our previous result using the RI/SMOM scheme [14,15], an improved version of the RI/MOM scheme with reduced sensitivity to long-distance effects, for the same ensembles [16].
An important future use for this sphere-averaged position-space renormalization scheme is to accurately define the weak operators which are needed in the calculation of non-leptonic decays, such as the K → ππ decay, in a three-flavor theory. At present these three-flavor operators are determined by using QCD perturbation theory to calculate that combination of three-flavor operators which will give the same matrix elements as the more physical four-flavor operators when evaluated at energies below the charm threshold. Such a use of QCD perturbation theory below the charm threshold introduces uncontrolled systematic errors. However, a nonperturbative matching of three-and four-flavor operators using RI/MOM methods is also potentially uncertain. The gauge-noninvariant operators that are traditionally neglected in RI/MOM calculations when performed at higher energies because of the presence of explicit factors of the gluon field, may give large contributions at energies below the charm mass. The sphere-averaged position-space renormalization scheme may allow a nonperturbative determination of the three-flavor Wilson coefficients in which the only errors, which are systematically improvable, come from the neglect of higher-dimension operators proportional to inverse powers of the charm quark mass.
The paper is organized as follows. In Section II, we summarize the traditional procedure of the position-space renormalization of an operator that does not mix with any other operator and identify the problem posed by the discretization errors that is addressed by the method presented in this paper. Our core technique in this work, the spherical average, is introduced in Section III. In Section IV, a concrete strategy to calculate the renormalized quark mass through the position-space renormalization of the scalar current is proposed.
In Section V, the details of the numerical simulation is described. In Section VI, our final result for the renormalized quark mass is shown. In the process, we show the performance of the spherical average especially at short distances and discuss how the renormalization window can be extended. In addition, we present a test of an ad hoc prescription to reduce nonperturbative effects at long distances that are mainly due to instanton interactions. In Section VII, we summarize the paper and discuss the prospect of further applications of the spherical average for various quantities calculated on the lattice. In Appendix A, we describe our investigation of the irregular a-dependence that appears in the spherical average, which turns out to be negligible.

II. FUNDAMENTAL PROCEDURE IN PREVIOUS WORKS
In this section, we summarize the traditional approach to position-space renormalization of an operator that does not mix with any other operator. We consider two-point Green's functions of a composite operator O s (µ; x) renormalized at a scale µ in a scheme s and the corresponding lattice operator O lat (1/a; an) for a lattice spacing a, Here, we distinguish a four-dimensional point in the continuum theory x = (x 1 , x 2 , x 3 , x 4 ) from that on the lattice an = (an 1 , an 2 , an 3 , an 4 ) since the discrete character of the lattice points is carefully considered throughout the paper. In this section, we treat these two-point functions in the chiral limit, which does not require consideration of the mass renormalization of quarks in the correlators and remove an extra scale from the renormalization procedure.
An operator O X (µ; x) renormalized at µ in the X-space scheme [2,10] is defined in the continuum theory by the condition where G free O (x) is the corresponding two-point Green's function evaluated in free field theory and |x| = µ x 2 µ . Since this nonperturbative scheme is fully gauge invariant and free from contact terms unlike the RI/MOM scheme, it prevents mixing with irrelevant operators and thus is a quite convenient scheme especially at low energies where perturbative schemes are not applicable and mixing with many irrelevant operators can occur in gauge noninvariant schemes.
The traditional renormalization condition yields the ensembles used to evaluate the continuum limit need to be integers or simple rational numbers. However, lattice spacings are not tuned so precisely in practical simulations. We propose a way to circumvent this problem in the next section.
We close the section by describing the relation between operators in the X-space scheme and those in another scheme s. Using Eqs. (3) and (4), the matching factor Z If we already know the correlator in the scheme s and any treatments in the X-space scheme such as the step scaling are not needed, we can skip renormalizing operators to the X-space scheme and directly compute Of course, this expression violates rotational symmetry as does Eq. (6) and therefore suffers from the same difficulty in taking the continuum limit of the corresponding renormalized operator as is described above.

III. SMOOTHING AVERAGE OVER SPHERES
Renormalization factors determined through the procedure discussed in the previous section contain discretization errors which depend in a complicated way on the lattice point n where the renormalization condition is imposed due to the violation of rotational symmetry. The complicated discretization errors induce difficulty in taking the continuum limit of renormalized operators as mentioned in the previous section. Some ideas to reduce this kind of discretization errors, subtracting free-field discretization error [2,11] and discarding the lattice data points where discretization errors are quite large [10,17], have been applied.
These previous works usually naïvely averaged the renormalization factor Eq.
with a coefficient c a,n which depends on the lattice point n in a complicated way. In the simplest case, F (x; a) is the continuum limit of the quantity being computed and does not depend on a. However, by including a possible logarithmic a-dependence, we can make our discussion more general and include the case where f a,n is an n-dependent renormalization factor such as the quantities given in Eqs. (6) and (8) or a correlator of unrenormalized operators.
We start with the case of one dimension, where we assume x = x 1 . We then use linear interpolation to extend the lattice results for f a,n , to define a functionf a (x) for all values of the continuous physical distance x: where n is now defined as x/a , the largest integer that is less than or equal to x/a.
Inserting Eq. (9) into this equation and expanding F (an; a) and F (a(n + 1); a) around x, we see thatf a (x) is an approximation to F (x; a) as a continuous function of x that is accurate up to O(a 2 ). Note that the appropriate weight of a(n + 1) − x and x − an in Eq. (10) is important to avoid introducing an O(a 1 ) error which would spoil the accuracy of an a 2 continuum extrapolation.
In the case of two dimensions, the weighted average Eq. (10) can be modified to a bilinear interpolation where n µ = x µ /a andμ is the unit vector for the µ-direction. While this weighted average is also easily found to be free from the O(a 1 ) error, it is expected to depend significantly on the direction of x as well as the distance |x| due to the violation of rotational symmetry.
The most naïve way to smooth this discretization error may be to introduce the average over a circle with the radius of |x|, where we use two-dimensional polar coordinates The extension to four dimensions is straightforward. The interpolation of f a at x is given byf where we define the factors One can easily verify this interpolated value is also free from the O(a 1 ) error. The smoothing average over the four-dimensional sphere with the radius of |x| iŝ with four-dimensional polar coordinates The averaged quantityf a (|x|) will differ from the direction independent continuum quantity Although the discretization error of the spherical average is thus O(a 2 ), it should be noted that the averaged value is not a regular polynomial in a but will contain extra nondifferentiable terms of O(a 2 ) because of the complicated a-dependence of the floor function n µ = x µ /a . This irregularity could arise also from the fact that the set of the lattice points n and their weight used by the spherical average at each fixed physical distance |x| depend on the lattice spacing a. Such complicated a-dependence could spoil the accuracy of a continuum extrapolation which assumed a regular a 2 term. In Appendix A, we discuss the significance of such complicated a-dependence and demonstrate it is small.

A. Strategy
Since the quark mass renormalization factor Z m can be calculated as the inverse of the renormalization factor Z S of the scalar current S(x) =ū(x)d(x), we consider the renormal- as long as chiral symmetry on the lattice is maintained. Since we use domain-wall fermions, we can calculate Z m from the renormalization of S(x) and P (x).
In what follows, we employ the MS scheme and introduce the input light quark mass parameter m ud that is used for the calculation of the correlators on the lattice. The ndependent renormalization factor Eq. (8) is then rewritten as 9 The chiral limit (m ud → 0) is taken in Section VI. In this work, the scalar correlator G MS S in continuum perturbation theory is considered only in the massless limit, where it is equivalent to the pseudoscalar correlator and is available to O((α s /π) 4 ) accuracy [12]. The strategy to improve the convergence of the perturbative series of the correlator is discussed in the following subsection and in [11] for more detail.
We also analyze an O(4)-symmetric renormalization factor obtained from Eq. (7) and the O(4)-symmetric renormalization condition with the sphere-averaged Green's function G lat S (1/a; |x|) calculated as follows. It should be noted that the complicated a-dependence appearing in the multi-linear interpolation depends on the first and second derivatives of the continuum version of the function that is to be interpolated with respect to |x|. Therefore, the spherical averaging procedure is applied to a function whose |x|-dependence in the continuum limit is as small as possible. For this reason, we calculate the spherical average of the ratio G lat S/P (1/a; an; m ud )/G free S (x)| x=an at each distance |x| and then define the sphere-averaged Green's function G lat S (1/a; |x|) as the product of it and G free S (x). Note that either Z MS/lat S/P (µ, 1/a; an; m ud ) or Z MS/lat S/P (µ, 1/a; |x|; m ud ) may not be an appropriate renormalization factor since it still depends on the location n or |x| due to the following sources of error: • Discretization effects in G lat S/P (1/a; an; m ud ).
• Truncation error from the perturbative calculation of G MS S (µ; x; 0).
• Nonperturbative QCD effects, which are not present in the perturbatively calculated G MS S (µ; x; 0) but do appear in the nonperturbatively measured G lat S/P (1/a; an; m ud ).
The first source is uncontrollable at short distances (|x|, a|n| ∼ a), while the others are significant at long distances (|x|, a|n| 1/Λ QCD ). We need to find or create an appropriate window where all of these sources of error are under control and the n-dependence of Using the unrenormalized quark mass m bare q (1/a) at the physical pion mass, which is given in Ref. [16] for the degenerate up and down quarks (q = ud) and the strange quark (q = s) on our ensembles, we analyze the n-and |x|-dependent renormalized quark masses and where q = ud, s.
In Section VI, we determine the renormalized mass of the degenerate up and down quarks and the strange quark on our ensembles.

B. Scalar correlator in massless perturbation theory
While the available four-loop perturbative results is an important advantage of the position-space renormalization of the scalar current, the region where discretization errors may be under controlled is 1/|x| 1 GeV for currently available lattices with domain-wall fermions and therefore the convergence of the perturbative expansion might be still insufficient. The convergence can be improved by a resummation of the perturbative series using the coupling constant at another renormalization scale as explained below.
Chetyrkin and Maier [12] gave the coefficients C S,CM i of the perturbative expansion up to i = 4. Here, the strong coupling constant a s (μ x ) = α s (μ x )/π is renormalized in the MS scheme atμ x = 2e −γ E /|x| 1.123/|x| with Euler's constant γ E = 0.5772 and is evaluated using the scale of QCD Λ MS QCD = 332(17) [18] in three flavor theory. By setting the renormalization scaleμ x of the scalar current and the strong coupling constant proportional to |x| −1 , the logarithmic |x|-dependence of the perturbative coefficients can be eliminated.
The anomalous dimension of the scalar current, which is the same as the mass anomalous dimension except for the sign and is calculated up to the five-loop level [19], enables us to evolve the scale on the LHS of Eq. (23). The beta function, which is also available to the five-loop level [20], can be used to evolve the scale of the strong coupling constant on the RHS of Eq. (23). Using the original perturbative coefficients C S,CM i and these scale evolution procedures, we obtain a general expression of the perturbative series where µ x and µ * x are the renormalization scale of the scalar current and that of the strong coupling constant, respectively. While the all-order calculation of the RHS is supposed to be independent of µ * x , any finite-order calculation does depend on µ * x . Therefore, the convergence of the perturbative series can be investigated by varying µ * x . Thus, we obtain the numerical value of the scalar correlator G MS with a scale µ * x of the strong coupling constant. In order to renormalize the scalar current at a specific scale, which we set to 3 GeV, the scale evolution is needed from µ x to 3 GeV, where γ m (z) and β(z) are the mass anomalous dimension and the beta function, respectively, and ρ(z) is known to the five-loop level [19]. While G MS S (3 GeV; x; 0) µ * x ,µ x is also supposed to be independent of µ * x and µ x in an all-order calculation, the convergence can be optimized by tuning the scale parameters µ x and µ * x so that the dependence on these scale parameters is minimized. We use the optimal values µ x = e 0.8 /|x| 2.2/|x| and µ * x = e 1.05 /|x| 2.9/|x| in the case of three-flavor QCD quoted by Ref. [11].

V. LATTICE SETUP
We perform lattice simulation with the ensembles of 2 + 1-flavor dynamical domain-wall fermions [21,22] and the Iwasaki gauge action [23,24] generated by the RBC and UKQCD collaborations [16]. Table I summarizes the properties of the ensembles used in this work.
We calculate with three lattice cutoffs ranging from 1.785(5) GeV to 3.148(17) GeV. The  (17) 1.005 (12) coarsest ensembles (24I) and the finer ones (32I and 32Ifine) are generated on the 24 3 × 64 and 32 3 × 64 lattices, respectively. The strange quark mass m s is used only for the strange sea quark, while we use the same values of the sea and valence quark masses m ud for the degenerate up and down quarks. The corresponding pion masses of the ensembles are quoted from Refs. [16,25] and in the region from 300 MeV to 430 MeV.
We distinguish the input quark masses m q (q = ud, s), which are used in the lattice calculations and shown in Table I, from the renormalized and unrenormalized quark masses (m MS q (µ) and m bare q (1/a)) that realize the physical pion mass. In Ref. [16], the values of unrenormalized quark masses were represented by where the values of quantities on the RHS were obtained by a global continuum and chiral fit to ten ensembles in the continuum scaling with the input experimental values of pion, kaon and Omega baryon masses [16]. We use m bare,32I  Table II.
and average the correlators over all these source points. The n-dependence of this quantity arises mainly from discretization errors at short distances, the truncation error of the perturbative calculation and nonperturbative effects at long distances as explained in Section IV A. These sources of the n-dependence need to be under controlled in order to obtain the correct value of the renormalized mass m MS q (3 GeV). However, Figure 1 indicates the ambiguity due to such n-dependence is O(1 MeV), which is much larger than the uncertainty of the renormalized light quark mass calculated by other works.
A rapid decrease is seen below 1/|x| = 1/a|n| ∼ 0.3 GeV since the truncation uncertainty of the perturbative calculation increases tremendously below this threshold. We do not expect that this lower limit on the perturbative window can be decreased because the convergence is already optimized by our choice of µ * x = e 1.05 /|x| and µ x = e 0.80 /|x|. These choices are found to maintain reasonable convergence down to 1/|x| ∼ 0.4 GeV [11].
Among the three sources of n-dependence listed in Section IV A, the n-dependence associated with the convergence of the perturbative calculation is thus already taken into account as much as possible. We discuss and take into account the remaining two sources below. The n-dependence associated with discretization errors can be reduced by the spherical average designed in Section III. The third source of n-dependence associated with nonperturbative effects can be investigated by comparing the scalar and pseudoscalar channels. While Figure 1 provides some information, we prefer to take the spherical average first and then to discuss the difference between the scalar and pseudoscalar correlators. cretization errors appear to be much reduced and the the n-dependence is made smaller, there must still be some dependence on the regularization parameter a in the form of a 2 /x 2 , a 2 Λ 2 QCD , a 4 /|x| 4 and so on. Therefore, the continuum extrapolation should be performed at sufficiently long distances where only O(a 2 ) discretization errors are visible. On the other hand, the difference between the scalar and pseudoscalar channels is significant, 1% at 1/|x| 1 GeV and 3% at 1/|x| 0.8 GeV, although we use domain-wall fermions which have high degree of chiral symmetry. Therefore, the extraction of the quark mass should be done at sufficiently short distances, 1/|x| 1 GeV in order to obtain a systematic error less than 1%. For this, the continuum extrapolation has to be safe at 1/|x| 1 GeV.
An important advantage of the spherical average is our ability to take the continuum limit of renormalized quantities as explained below. Since the structure of the a-dependence depends on |x| as it contains a term proportional to a 2 /x 2 , the renormalized quantities need to be calculated at the same physical distance scale |x| for each ensemble in order with three fit parameters: m MS ud,S/P (3 GeV; |x|), C a,S/P (|x|) and C m,S/P (|x|) for each |x|. Here, we introduce a term proportional to the pion mass squared M π (a, m ud ) 2 labeled by the ensemble parameters a and m ud , although the leading mass correction in perturbation theory is proportional to quark mass squared or M π (a, m ud ) 4 . This is because the perturbative mass correction is much smaller than the mass correction from OPE such as m qq and m qGq around 1/|x| ∼ 0.5 GeV [17,26].    Table I and its extrapolation to the continuum and chiral limits (a → 0, m ud → 0).
The results for the scalar (upper panel) and pseudoscalar (lower panel) channels are shown separately.
ensemble listed in Table I    We show the result only for 1/|x| ≤ 0.9 GeV because the spherical average with |x| < 2a, which corresponds to 1/|x| > 0.89 GeV on the coarsest lattice, uses the value at x = 0 and therefore contains unphysical contact terms. Since there is no plateau in the extrapolated results and the difference between the scalar and pseudoscalar is quite significant, it is difficult to determine the quark mass from these plots. We may need to exclude the data on the coarsest lattice from this analysis. Figure 4 shows the results on the finer two lattices (32I and 32Ifine). The continuum limit is taken only with these lattice data, excluding the result on the coarsest ensembles. Since the number of free parameters in Eq. (28) and the number of data points in this extrapolation are both three, the extrapolation is not an actual χ 2 fit but the free parameters can still be determined, with propagated errors, by solving Eq. (28). While the |x|-dependence of the extrapolated result becomes milder especially around 1/|x| 0.8 GeV, the statistical error is substantially increased by discarding the data on the coarsest lattice. In order to obtain a reasonable result with sufficiently small statistical error from such an analysis, we need to introduce finer lattices.
Since it is currently not easy to introduce a finer lattice, we seek a more economical analysis that enables to extract the quark mass from the data we currently have. Among the three sources of the n-dependence of Z MS/lat S/P (µ, 1/a; an; m ud ) mentioned in Section IV A, we now focus on the third one, the nonperturbative effects. The nonperturbative effects on the scalar and pseudoscalar correlators are known to be quite large compared to those on the vector and axial-vector correlators in a model based on instantons because the scalar and pseudoscalar channels are directly affected by instantons [27,28]. The effect of a single instanton on these channels, which is the most significant at short distances, is of the same magnitude but with opposite sign [29]. Therefore, the naïve average of these two channels may be free from the largest source of nonperturbative effects. Therefore, we analyze with the O(4)-symmetric renormalization factor Z MS/lat S+P (3 GeV, 1/a; |x|; m ud ) obtained from Eqs. (19) and (20) by substituting S/P with S + P and defining G lat S+P (1/a; |x|) as the spherical average of the average of the scalar and pseudoscalar correlators.
with the fit parameters m        To investigate the a-and M π -dependences of m MS ud (3 GeV; |x|; a, m ud ) more clearly, we visualize this extrapolation by analyzing the renormalized mass at a specific distance |x| in the chiral limit and in the continuum limit with the parameters m MS ud (3 GeV; |x|), C a (|x|) and C m (|x|) obtained through the simultaneous fit Eq. (30).
The first error is the statistical error. The second error is the systematic error due to discretization effects, which is estimated by increasing 1/|x| up to 0.60 GeV where a deviation from the plateau is beginning. The third error is the systematic uncertainty due to the truncation of the perturbative calculation, which is estimated by varying the parameter µ * x in the region µ * ,opt  [30] methods using 2 + 1-flavor ensembles. Applying the same procedure, we obtain the strange  [16] determined through the RI/SMOM scheme using the same ensembles.
In Table III, we summarize the quark mass renormalization factors calculated on each ensemble and at three values of |x|. The three errors shown (from left to right) correspond to the statistical error, the systematic error due to the truncation of the perturbative calculation and due to the uncertainty of the QCD scale Λ MS QCD , in order. The latter two uncertainties, which are associated with perturbation theory, are more significant at longer distances as one can easily expect. As we have mentioned throughout the paper, we need to fix |x| when we renormalize a quantity and take its continuum and chiral limits since the structure of aand m ud -dependences depends on |x|.

VII. CONCLUSION
We have proposed a spherical averaging technique for position-space renormalization to reduce discretization errors and enhance the renormalization window. This technique has the further important advantage that it allows renormalized quantities to be defined at any fixed physical distance. This allows a direct matching between renormalized quantities defined on ensembles with different lattice spacings and a continuum limit to be easily taken for position-space renormalized quantities at a fixed, physical renormalization scale.
The technique is applied to the quark mass renormalization using the scalar and pseudoscalar correlators in position space. We investigate the |x|-dependence of the renormalized quark mass in the continuum limit and find a plateau even when the MS renormalized quark mass at finite a still depends slightly on the distance |x| at which the intermediate position- imposed. These will be insufficient to determine the needed N × N renormalization matrix.
An extension of the spherical averaging procedure to such three-point functions is the next step with is being developed.  9. B 1 (x; a)/(x 2 F (x)), B 2 (|x|; a)/(x 2 F rr (|x|)) and B 4 (|x|; a)/(x 2 F rr (|x|)) plotted as functions of a/|x|. The curve of (a/|x|) 2 /8 is also plotted for comparison.
While the both of the second and third terms in Eq. (A1) are expected to contain complicated a-dependence, the third term can be explicitly analyzed and therefore is discussed first. Since the value of n = x/a jumps where x/a is an integer, the a-dependence of I 1 (x; a)/(x 2 F (x)) given by Eq. (A2) is drawn (dashed curve) in Figure 9. Thus, the continuum extrapolation with a few data points with an assumption of simple a 2 discretization error could be inaccurate. While such ambiguity is expected to be less than 1% of x 2 F (x) in the case of one dimension, one could imagine that in higher dimensions the spherical average further softens such irregular dependence on a and makes the continuum extrapolation more accurate.
The interpolation in d dimensions can be written as where and F µµ (x) is the second derivative of F (x) with respect to x µ . Averaging over the sphere by the integral given in Eq. (12) 10. B 2 (|x|; a)/(a 2 F rr (|x|)) and B 4 (|x|; a)/(a 2 F rr (r)) shifted by δ d , the mid point of the oscillation. Note by dividing by a 2 instead of by |x| 2 as is done in Figure 9, we are plotting the correction relative to the regular a 2 error.
term becomes where and F r (|x|) and F rr (|x|) respectively stand for the first and second derivatives of F (|x|) with respect to |x|.
In order to describe the a dependence implied by Eq. (A5), we need to assume a relation between the F rr (|x|) and F r (|x|)/|x| terms found in the first line of that equation. For the purposes of illustration we will assume that F (|x|), which in this application is expected to be a slowly varying function of |x|, behaves as ln(|x|) so we can replace F r (|x|) by the |x|F rr (|x|). In this case, Eq. (A5) becomes  smaller than that in one dimension and therefore the continuum extrapolation with a regular a 2 term is expected to be more accurate. The irregular a-dependence of the spherical average in four dimensions is even smaller than that in two dimensions in the sense explained below.
The significance of the irregular term in B d (|x|; a) x 2 F rr (|x|) = δ d · (a/|x|) 2 + (irregular oscillation) + O (a/|x|) 3 , can be investigated by which is plotted in Figure 10 for two and four dimensions. Here, we find δ 2 0.00047 and δ 4 0.16667. As the figure indicates, the irregular oscillation in four dimensions is even smaller than that in two dimensions. It means the continuum limit can be safely taken with a regular a 2 term.
Thus, we conclude that the irregular a-dependence associated with the third term of Eq. (A1) or (A3) is negligible. We proceed to discuss the other source of the irregular a-dependence associated with c a (x) in Eq. (A1) or (A3), which can be written as c a (x) = a −4 1 i,j,k,l=0 ∆ 1,i ∆ 2,j ∆ 3,k ∆ 4,l c a,n+i1+j2+k3+l4 , in the case of four dimensions. Here, c n for the scalar or pseudoscalar correlator may be roughly approximated to a dispersion integral of the difference between the lattice and continuum propagators of a scalar field c a,n a 2 = The result is shown in Figure 11, which indicates that the discretization error is mostly proportional to a 2 for |x| 3a and that the continuum extrapolation using lattice data at |x| 3a, 4a and 5a is likely accurate within the O(0.1%) level.
While the above analysis using the bosonic propagator may be valid in QCD at long distances, the discretization error of Green's functions in the perturbative regime may need to be discussed in terms of fermionic propagators. We analyze the spherical average∆ 2f (|x|; 0) of ∆ 2f (an; 0) = lim as the Fourier transform of the corresponding momentum-space propagator of domain-wall fermions [22,31]. The calculation on a 200 4 lattice at am ∼ 0.005 still suffers from significant finite volume effects, which can mostly be estimated as the effect of the fermions wrapping around the volume in the continuum theory. See [11] for more detail. Figure 12 shows the result for the spherical average∆ 2f (|x|; 0) plotted as a function of (a/|x|) 2 (upper panel) and (a/|x|) 4 (lower panel).Although the a-dependence at small values of a appears to be proportional to a 4 , the coefficient is somehow large ∼ 50(a/|x|) 4 . We might therefore need to be careful when we take the continuum limit. Of course, nonperturbative interactions