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

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

Masaaki Tomii    Norman H. Christ Physics Department, Columbia University, New York 10027, USA
###### Abstract

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 -symmetric position-space Green’s functions in the continuum theory but with a residual dependence on a regularization parameter, the lattice spacing . We can then take the continuum limit of these renormalized quantities calculated at the same physical renormalization scale and investigate the resulting -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 -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  Martinelli:1997zc (); Gimenez:2004me (), 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 Martinelli:1994ty (). 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  Luscher:1991wu (); Jansen:1995ck (); Arthur:2010ht (). (For example, when Cichy et al. Cichy:2016qpu () 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 Symanzik:1983dc (); Symanzik:1983gh (). 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

 G(xn)=⟨O(xn)O(0)†⟩ (1)

where is a gauge-invariant composite local operator, is a point on our discrete lattice determined by the four integers and, for simplicity, the second point 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 of the continuous position four-vector , obtained by multi-linear interpolation from the sixteen values obtained by evaluating at the sixteen lattice points that lie at the vertices of the four-dimensional cube which contains the point .

Assuming, as we do throughout this paper, that our lattice theory has no order errors, our interpolated Green’s function will agree with the corresponding Green’s function of the continuum theory up to errors which vanish as in the continuum limit. Of course, the errors which appear in the piece-wise linear function will still reflect the 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 over a three-dimension sphere of radius centered at the origin:

 ˆG(|x|)=12π2|x|3∫d4x′δ(|x′|−|x|)¯¯¯¯G(x′). (2)

We will then impose conditions on to renormalize the operator .

The finite lattice spacing errors that are present in the lattice Green’s function are expected to appear as simple polynomials in with coefficients which in perturbation theory depend only logarithmically in , allowing a simple exptrapolation to the continuum limit. The lattice spacing dependence of our averaged quantity will be more complicated. In addition to the simple errors coming from , the sphere averaging procedure will introduce errors which, while bounded by may be complicated irregular functions of which could cause an explicit extrapolation in to fail. As will be shown in Appendix A, this irregular dependence on 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 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  Gimenez:2004me (); Cichy:2012is (); Tomii:2016xiv (). For renormalization of bilinear operators, there is another important advantage of the position-space procedure: the perturbative matching to the modified minimal subtraction () scheme is available to for the vector, axial-vector, scalar and pseudoscalar currents and to for the tensor current  Chetyrkin:2010dx (). Utilizing the spherical averaging technique, we perform a new analysis that takes the continuum limit of the renormalized quark mass at many values of and shows its -dependence. The final result agrees with the FLAG average Aoki:2016frl () as well as our previous result using the RI/SMOM scheme  Aoki:2007xm (); Sturm:2009kb (), an improved version of the RI/MOM scheme with reduced sensitivity to long-distance effects, for the same ensembles Blum:2014tka ().

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 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 -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 renormalized at a scale in a scheme and the corresponding lattice operator for a lattice spacing ,

 GsO(μ;x)=⟨Os(μ;x)Os(μ;0)†⟩,  GlatO(1/a;an)=⟨Olat(1/a;an)Olat(1/a;0)†⟩. (3)

Here, we distinguish a four-dimensional point in the continuum theory from that on the lattice 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 renormalized at in the X-space scheme Gimenez:2004me (); Cichy:2012is () is defined in the continuum theory by the condition

 GXO(μ;x)∣∣μ=1/|x|=GfreeO(x), (4)

where is the corresponding two-point Green’s function evaluated in free field theory and . 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.

 ˜ZX/latO(μ,1/a;an)2∣∣μ=1/a|n|GlatO(1/a;an)=GfreeO(x)∣∣x=an, (5)

yields

 ˜ZX/latO(μ,1/a;an)∣∣μ=1/a|n|= ⎷GfreeO(x)∣∣x=anGlatO(1/a;an), (6)

which violates rotational symmetry and depends on in a complicated way. Since the -violating -dependence is , it can be eliminated and only the dependence on the distance scale remains if the continuum limit of the renormalized operator is accurately taken. However, evaluating the continuum limit requires an extrapolation of numerical values at a fixed physical location so that when comparing ensembles and it is only the lattice spacing, not the physical position which is changing. This means the ratios of the lattice spacings for 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 -space scheme and those in another scheme . Using Eqs. (3) and (4), the matching factor , which is defined by , can be written as

 Zs/XO(μ,μ′)= ⎷GsO(μ;x)GfreeO(x)∣∣∣|x|=1/μ′. (7)

If we already know the correlator in the scheme 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

 ˜Zs/latO(μ,1/a;an)≡Zs/XO(μ,μ′)˜ZX/latO(μ′,1/a;an)∣∣μ′=1/a|n|= ⎷GsO(μ;x)∣∣x=anGlatO(1/a;an). (8)

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 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 Gimenez:2004me (); Tomii:2016xiv () and discarding the lattice data points where discretization errors are quite large  Chu:1993cn (); Cichy:2012is (), have been applied. These previous works usually naïvely averaged the renormalization factor Eq. (8) over lattice points in the renormalization window, which could induce an irrelevant linear dependence on and further degrade the accuracy of the continuum extrapolation of a renormalized quantity which assumed that the leading discretization error is . In this section, we propose another way to smooth lattice results, in which the irrelevant discretization error does not appear and the continuum extrapolation of a renormalized quantity using a constant plus an term can be safely taken.

We consider a lattice quantity calculated at each lattice point . The -dependence of can be sketched as

 fa,n=F(x;a)|x=an+ca,na2+O(a4), (9)

with a coefficient which depends on the lattice point in a complicated way. In the simplest case, is the continuum limit of the quantity being computed and does not depend on . However, by including a possible logarithmic -dependence, we can make our discussion more general and include the case where is an -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 . We then use linear interpolation to extend the lattice results for , to define a function for all values of the continuous physical distance :

 ¯fa(x)=(a(n+1)−x)fa,n+(x−an)fa,n+1a, (10)

where is now defined as , the largest integer that is less than or equal to . Inserting Eq. (9) into this equation and expanding and around , we see that is an approximation to as a continuous function of that is accurate up to . Note that the appropriate weight of and in Eq. (10) is important to avoid introducing an error which would spoil the accuracy of an continuum extrapolation.

In the case of two dimensions, the weighted average Eq. (10) can be modified to a bilinear interpolation

 ¯fa(x) (11)

where and is the unit vector for the -direction. While this weighted average is also easily found to be free from the error, it is expected to depend significantly on the direction of as well as the distance 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 ,

 ^fa(|x|)=12π∫2π0dθ¯fa(x), (12)

where we use two-dimensional polar coordinates

 x1=|x|cosθ,  x2=|x|sinθ. (13)

The extension to four dimensions is straightforward. The interpolation of at is given by

 ¯fa(x)=a−41∑i,j,k,l=0Δ1,iΔ2,jΔ3,kΔ4,lfa,n+i^1+j^2+k^3+l^4, (14)

where we define the factors

 Δμ,i=|a(nμ+1−i)−xμ|. (15)

One can easily verify this interpolated value is also free from the error. The smoothing average over the four-dimensional sphere with the radius of is

 ^fa(|x|)=12π2∫π0dθ1∫π0dθ2∫2π0dθ3sin2θ1sinθ2¯fa(x), (16)

with four-dimensional polar coordinates

 x1 =|x|cosθ1, x2 =|x|sinθ1cosθ2, x3 =|x|sinθ1sinθ2cosθ3, x4 =|x|sinθ1sinθ2sinθ3. (17)

The averaged quantity will differ from the direction independent continuum quantity by discretization errors of .

Although the discretization error of the spherical average is thus , it should be noted that the averaged value is not a regular polynomial in but will contain extra non-differentiable terms of because of the complicated -dependence of the floor function . This irregularity could arise also from the fact that the set of the lattice points and their weight used by the spherical average at each fixed physical distance depend on the lattice spacing . Such complicated -dependence could spoil the accuracy of a continuum extrapolation which assumed a regular term. In Appendix A, we discuss the significance of such complicated -dependence and demonstrate it is small.

## Iv Quark masses renormalization in position space

### iv.1 Strategy

Since the quark mass renormalization factor can be calculated as the inverse of the renormalization factor of the scalar current , we consider the renormalization of , which is equivalent to that of the pseudoscalar current as long as chiral symmetry on the lattice is maintained. Since we use domain-wall fermions, we can calculate from the renormalization of and .

In what follows, we employ the scheme and introduce the input light quark mass parameter that is used for the calculation of the correlators on the lattice. The -dependent renormalization factor Eq. (8) is then rewritten as

 ˜Z¯¯¯¯¯¯¯MS/latS/P(μ,1/a;an;m′ud)=  ⎷G¯¯¯¯¯¯¯MSS(μ;x;0)∣∣x=anGlatS/P(1/a;an;m′ud). (18)

The chiral limit () is taken in Section VI. In this work, the scalar correlator in continuum perturbation theory is considered only in the massless limit, where it is equivalent to the pseudoscalar correlator and is available to accuracy Chetyrkin:2010dx (). The strategy to improve the convergence of the perturbative series of the correlator is discussed in the following subsection and in Tomii:2016xiv () for more detail.

We also analyze an -symmetric renormalization factor

 ˆ˜Z¯¯¯¯¯¯¯MS/latS/P(μ,1/a;|x|;m′ud)=Z¯¯¯¯¯¯¯MS/XS(μ,μ′)ˆ˜ZX/latS/P(μ′,1/a;|x|;m′ud) (19)

obtained from Eq. (7) and the -symmetric renormalization condition

 (20)

with the sphere-averaged Green’s function calculated as follows. It should be noted that the complicated -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 . Therefore, the spherical averaging procedure is applied to a function whose -dependence in the continuum limit is as small as possible. For this reason, we calculate the spherical average of the ratio at each distance and then define the sphere-averaged Green’s function as the product of it and .

Note that either or may not be an appropriate renormalization factor since it still depends on the location or due to the following sources of error:

• Discretization effects in .

• Truncation error from the perturbative calculation of .

• Nonperturbative QCD effects, which are not present in the perturbatively calculated but do appear in the nonperturbatively measured .

The first source is uncontrollable at short distances (), while the others are significant at long distances (). We need to find or create an appropriate window where all of these sources of error are under control and the -dependence of or -dependence of is sufficiently small. Since the third source especially violates the degeneracy of and , analyzing both of these may specify the region where nonperturbative effects are less significant.

Using the unrenormalized quark mass at the physical pion mass, which is given in Ref. Blum:2014tka () for the degenerate up and down quarks and the strange quark on our ensembles, we analyze the - and -dependent renormalized quark masses

 ˜m¯¯¯¯¯¯¯MSq,S/P(μ;an;a,m′ud)=mbareq(1/a)˜Z¯¯¯¯¯¯¯MS/latS/P(μ,1/a;an;m′ud), (21)

and

 ˆ˜m¯¯¯¯¯¯¯MSq,S/P(μ;|x|;a,m′ud)=mbareq(1/a)ˆ˜Z¯¯¯¯¯¯¯MS/latS/P(μ,1/a;|x|;m′ud), (22)

where .

In Section VI, we determine the renormalized mass of the degenerate up and down quarks and the strange quark on our ensembles.

### iv.2 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  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 Chetyrkin:2010dx () gave the coefficients of the perturbative expansion

 G¯¯¯¯¯¯¯MSS(~μx;x;0)=3π4|x|6(1+∑iCS,CMias(~μx)i), (23)

up to . Here, the strong coupling constant is renormalized in the scheme at with Euler’s constant and is evaluated using the scale of QCD  Tanabashi:2018oca () in three flavor theory. By setting the renormalization scale of the scalar current and the strong coupling constant proportional to , the logarithmic -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 Baikov:2014qja (), enables us to evolve the scale on the LHS of Eq. (23). The beta function, which is also available to the five-loop level Baikov:2016tgj (), can be used to evolve the scale of the strong coupling constant on the RHS of Eq. (23). Using the original perturbative coefficients and these scale evolution procedures, we obtain a general expression of the perturbative series

 G¯¯¯¯¯¯¯MSS(μ′x;x;0)=3π4x6(1+∑iCSi(μ∗x,μ′x)as(μ∗x)i), (24)

where and 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 , any finite-order calculation does depend on . Therefore, the convergence of the perturbative series can be investigated by varying .

Thus, we obtain the numerical value of the scalar correlator calculated with a scale 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 to 3 GeV,

 G¯¯¯¯¯¯¯MSS(3GeV;x;0)∣∣μ∗x,μ′x =exp(−2∫as(3GeV)as(μ′x)dzzγm(z)β(z))G¯¯¯¯¯¯¯MSS(μ′x;x;0)∣∣μ∗x =(ρ(as(μ′x))ρ(as(3GeV)))2G¯¯¯¯¯¯¯MSS(μ′x;x;0)∣∣μ∗x, (25)

where and are the mass anomalous dimension and the beta function, respectively, and is known to the five-loop level  Baikov:2014qja (). While is also supposed to be independent of and in an all-order calculation, the convergence can be optimized by tuning the scale parameters and so that the dependence on these scale parameters is minimized. We use the optimal values and in the case of three-flavor QCD quoted by Ref. Tomii:2016xiv ().

## V Lattice setup

We perform lattice simulation with the ensembles of -flavor dynamical domain-wall fermions Kaplan:1992bt (); Shamir:1993zy () and the Iwasaki gauge action Iwasaki:1984cj (); Iwasaki:1985we () generated by the RBC and UKQCD collaborations Blum:2014tka (). Table 1 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 coarsest ensembles (24I) and the finer ones (32I and 32Ifine) are generated on the and lattices, respectively. The strange quark mass is used only for the strange sea quark, while we use the same values of the sea and valence quark masses for the degenerate up and down quarks. The corresponding pion masses of the ensembles are quoted from Refs. Aoki:2010dy (); Blum:2014tka () and in the region from 300 MeV to 430 MeV.

We distinguish the input quark masses , which are used in the lattice calculations and shown in Table 1, from the renormalized and unrenormalized quark masses ( and ) that realize the physical pion mass. In Ref. Blum:2014tka (), the values of unrenormalized quark masses were represented by

 mbareq(1/a)=mbare,32IqZq(1/a), (26)

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 Blum:2014tka (). We use  MeV,  MeV and summarized in Table 2.

For each configuration, we calculate the scalar and pseudoscalar correlators with 16 point sources located at

 (0,0,0,0),  (L2,L2,0,0),  (0,L2,L2,0),  (L2,0,L2,0),(L4,L4,L4,T4),  (3L4,3L4,L4,T4),  (L4,3L4,3L4,T4),  (3L4,L4,3L4,T4),(L2,L2,L2,T2),  (0,0,L2,T2),  (L2,0,0,T2),  (0,L2,0,T2),(3L4,3L4,3L4,3T4),  (L4,L4,3L4,3T4),  (3L4,L4,L4,3T4),  (L4,3L4,L4,3T4), (27)

and average the correlators over all these source points.

## Vi Numerical results

Figure 1 shows defined in Eq. (21) calculated on the 32Ifine ensemble. The results for both the scalar (diamonds) and pseudoscalar (circles) channels are shown. Because of the violation of rotational symmetry, we distinguish, in the figure, different lattice points that are not equivalent with respect to rotations or parity inversion in the four-dimensional hypercubic group. For example, (1,1,1,1) and (2,0,0,0) correspond to the same distance but are distinguished since they are indeed different points if rotational symmetry is violated and only hypercubic symmetry remains. The results are averaged over sets of lattice points related by rotations including parity inversion.

The -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.1. These sources of the -dependence need to be under controlled in order to obtain the correct value of the renormalized mass . However, Figure 1 indicates the ambiguity due to such -dependence is , which is much larger than the uncertainty of the renormalized light quark mass calculated by other works.

A rapid decrease is seen below  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 and . These choices are found to maintain reasonable convergence down to  GeV Tomii:2016xiv ().

Among the three sources of -dependence listed in Section IV.1, the -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 -dependence associated with discretization errors can be reduced by the spherical average designed in Section III. The third source of -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.

Figure 2 shows the result for defined in Eq. (22), where the renormalization factors of the scalar (diamonds) and pseudoscalar (circles) currents are calculated using the spherical average of the corresponding Green’s functions calculated on the 32I ensembles with (filled points) and (open points). While discretization errors appear to be much reduced and the the -dependence is made smaller, there must still be some dependence on the regularization parameter in the form of , , and so on. Therefore, the continuum extrapolation should be performed at sufficiently long distances where only discretization errors are visible. On the other hand, the difference between the scalar and pseudoscalar channels is significant, 1% at  GeV and 3% at  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,  GeV in order to obtain a systematic error less than 1%. For this, the continuum extrapolation has to be safe at  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 -dependence depends on as it contains a term proportional to , the renormalized quantities need to be calculated at the same physical distance scale for each ensemble in order to take the continuum limit as a quadratic extrapolation. The spherical averaging technique trivially enables such an extrapolation. The extrapolation of the spherical average to the continuum () and chiral () limits is done by performing a simultaneous fit to the data from all the ensembles with the fit function

 ˆ˜m¯¯¯¯¯¯¯MSud,S/P(3 GeV;|x|;a,m′ud)=ˆ˜m¯¯¯¯¯¯¯MSud,S/P(3 GeV;|x|)+Ca,S/P(|x|)a2+Cm,S/P(|x|)Mπ(a,m′ud)2, (28)

with three fit parameters: , and for each . Here, we introduce a term proportional to the pion mass squared labeled by the ensemble parameters and , although the leading mass correction in perturbation theory is proportional to quark mass squared or . This is because the perturbative mass correction is much smaller than the mass correction from OPE such as and around  GeV Chu:1993cn (); Hands:1994cj ().

Figure 3 shows the result for the spherical average on each ensemble listed in Table 1 and its continuum and chiral limits . We also show the FLAG average Aoki:2016frl () (solid band) and our previous RBC/UKQCD result obtained through the RI/SMOM scheme Blum:2014tka () (hatched band), which is  MeV including the statistical and systematic errors. While FLAG gave the value renormalized at 2 GeV Aoki:2016frl (), we perform its scale evolution to 3 GeV and show  MeV in the figure. We show the result only for  GeV because the spherical average with , which corresponds to  GeV on the coarsest lattice, uses the value at 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 fit but the free parameters can still be determined, with propagated errors, by solving Eq. (28). While the -dependence of the extrapolated result becomes milder especially around  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 -dependence of mentioned in Section IV.1, 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 'tHooft:1976fv (); Novikov:1981xi (). 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 Shuryak:1993kg (). Therefore, the naïve average of these two channels may be free from the largest source of nonperturbative effects. Therefore, we analyze

 ˆ˜m¯¯¯¯¯¯¯MSq(3 GeV;|x|;a,m′ud)=mbareq(1/a)ˆ˜Z¯¯¯¯¯¯¯MS/latS+P(3 GeV,1/a;|x|;m′ud), (29)

with the -symmetric renormalization factor obtained from Eqs. (19) and (20) by substituting with and defining as the spherical average of the average of the scalar and pseudoscalar correlators.

Figure 5 shows the results for . The continuum and chiral limits are taken using the fit function

 ˆ˜m¯¯¯¯¯¯¯MSud(3 GeV;|x|;a,m′ud)=ˆ˜m¯¯¯¯¯¯¯MSud(3 GeV;|x|)+Ca(|x|)a2+Cm(|x|)Mπ(a,m′ud)2, (30)

with the fit parameters , and . The fit results with (upper panel) and without (lower panel) the data on the coarsest lattice are shown. We see a plateau of the extrapolated data in the interval  GeV in the upper panel and  GeV in the lower panel. These facts agree with the instanton-based observation that the nonperturbative effects on the average of the scalar and pseudoscalar correlators are much smaller than those on the individual channels.

Figure 6 shows the result for the strange quark mass defined in Eq. (29) with , where the same renormalization factors as for the light quark mass are used. The continuum and chiral extrapolations are done by the fit function (30) with the substitution . A plateau is seen in the same region as in the result for the light quarks mass.

Figure 7 shows the d.o.f. obtained through the simultaneous fit both for the light (crosses) and strange (circles) quark masses. While the position-space renormalization factors calculated on each ensemble are uncorrelated, the statistical errors for and , which are taken from Ref. Blum:2014tka () and used in Eq. (26) to calculate unrenormalized quark mass for each lattice cutoff, are likely correlated. We interpret the small values of d.o.f. shown in Figure 6 as resulting from our ignorance of such correlations. Since /d.o.f. at  GeV, which roughly corresponds to on the coarsest ensemble, is not too large, we may conclude that the spherical average does not suffer significantly from higher orders of errors for .

To investigate the - and -dependences of more clearly, we visualize this extrapolation by analyzing the renormalized mass at a specific distance in the chiral limit

 ˆ˜m¯¯¯¯¯¯¯MSud(3 GeV;|x|;a,m′ud→0) ≡ˆ˜m¯¯¯¯¯¯¯MSud(3 GeV;|x|;a,m′ud)−CmMπ(a,m′ud)2 =