Perturbative and Nonperturbative Renormalization in Lattice QCD
We investigate the perturbative and nonperturbative renormalization of composite operators in lattice QCD restricting ourselves to operators that are bilinear in the quark fields (quark-antiquark operators). These include operators which are relevant to the calculation of moments of hadronic structure functions. The nonperturbative computations are based on Monte Carlo simulations with two flavors of clover fermions and utilize the Rome-Southampton method also known as the RI-MOM scheme. We compare the results of this approach with various estimates from lattice perturbation theory, in particular with recent two-loop calculations.
The investigation of hadron structure has become a central topic of lattice QCD. In many cases this involves the evaluation of matrix elements of local operators between hadron states. For example, (moments of) generalized parton distributions can be extracted from matrix elements of quark-antiquark operators, i.e., operators composed of a quark field, its adjoint and a number of gluon fields entering through covariant derivatives which act on the quark fields. In general such operators have to be renormalized. In this process the operator of interest may receive contributions also from other operators, i.e., it may mix with these additional operators. On the lattice, mixing occurs more frequently than in the continuum due to the restricted space-time symmetries. Since in the end one wants to make contact with phenomenological studies, which almost exclusively refer to operators renormalized in the scheme of dimensional regularization, one needs the renormalization factors leading from the bare operators on the lattice to the operators in the continuum.
The most straightforward approach towards the calculation of renormalization factors is based on lattice perturbation theory (for a review see Ref. stefano ()). Unfortunately, this method meets with some difficulties. First, perturbation theory on the lattice is computationally much more complex than in the continuum and therefore the calculations rarely extend beyond one-loop order (see, however, Refs. mason (); panagopoulos1 (); panagopoulos2 ()). Second, lattice perturbation theory usually converges rather slowly so that the accuracy of perturbative renormalization factors is limited. Identifying one source of these poor convergence properties, Lepage and Mackenzie proposed as a remedy the so-called tadpole improved perturbation theory lepmac (). Still, considerable uncertainty remains. Third, mixing with operators of lower dimension cannot be treated by perturbation theory.
In special cases, when the renormalization factors contain no ultraviolet divergences, a nonperturbative determination is possible with the help of Ward identities bochi (). However, there are many interesting operators that cannot be renormalized by this method.
A general nonperturbative approach to renormalization has been developed within the Schrödinger functional scheme (see, e.g., Refs. sf1 (); sf2 (), reviews are given in Refs. sommer1 (); sommer2 ()). In this method the finite size of the lattices employed in the simulations (along with appropriate boundary conditions in Euclidean time) is used to set the renormalization scale. In the end continuum perturbation theory is employed to convert the results from the Schrödinger functional scheme to the scheme. Though theoretically appealing the practical implementation of the procedure requires a lot of effort and has to be repeated for every new operator again from the very beginning.
Another nonperturbative method for computing renormalization coefficients of arbitrary quark-antiquark operators is the Rome-Southampton method (also known as the RI-MOM scheme) introduced in Ref. marti (). It mimics the procedure used in continuum perturbation theory. The basic objects are quark two-point functions with an insertion of the operator under consideration at momentum zero. These are computed in a suitable gauge, e.g., the Landau gauge. In continuum perturbation theory the two-point functions are calculated order by order in an expansion in powers of the strong coupling constant while in the Rome-Southampton method they are evaluated within a Monte Carlo simulation on the lattice. In order to extract from these data renormalization factors which yield renormalized operators in the scheme in the continuum limit one needs a renormalization condition which is applicable to lattice data as well as to perturbative continuum results. A suitable condition has been given in Ref. marti ().
Compared with the Schrödinger functional approach the Rome-Southampton method is distinguished by its relatively simple implementation. Furthermore, one can deal with all desired operators in a single simulation. On the other hand, the Schrödinger functional method is explicitly gauge invariant, while the Rome-Southampton method requires gauge fixing.
In a previous publication reno () we have performed an extensive study of nonperturbative renormalization for a variety of quark-antiquark operators using the Rome-Southampton method, motivated by our investigations of hadron structure functions. This was done with Wilson fermions in the quenched approximation for two values of the lattice spacing. Later on, these studies have been extended to improved Wilson fermions, based on quenched simulations at three values of the lattice spacing timid (). Meanwhile we are using gauge field configurations generated with two flavors of dynamical quarks, which has made a reconsideration of renormalization necessary.
In this paper we present results for renormalization factors obtained within the Rome-Southampton approach with dynamical quarks. We work with nonperturbatively -improved Wilson fermions (clover fermions). The operators are, however, not (yet) improved. We continue to apply the momentum sources introduced in Ref. reno (). In addition, we have refined our approach, subtracting lattice artefacts through one-loop boosted perturbation theory. As in Ref. reno () we consider only flavor-nonsinglet quark-antiquark operators. For some thoughts concerning flavor-singlet operators see Refs. qmass (); fermilab ().
The paper is organized as follows: After introducing in Sec. II the operators to be studied we explain the method of nonperturbative renormalization in Sec. III. Our implementation of this method employing momentum sources is described in Sec. IV. After a brief overview over our gauge field configurations in Sec. V we discuss the chiral extrapolation of our data in Sec. VI. Section VII reviews formulae from continuum perturbation theory that will be needed in the analysis. Results from lattice perturbation theory are compiled in Sec. VIII. Section IX explains how we apply lattice perturbation theory in order to subtract lattice artefacts. In Sec. X we describe our method of extracting the renormalization factors from the Monte Carlo data. The results (perturbative as well as nonperturbative) are then presented and discussed in Sec. XI. Finally, we present our conclusions in Sec. XII. Some technical details are explained in the Appendices.
Ii The operators
In the Euclidean continuum we want to study the operators
(with and ) or rather O(4) irreducible multiplets with definite charge conjugation parity. In particular, we obtain twist-2 operators by symmetrizing the indices and subtracting the traces. We have given the quark fields definite flavors (assumed to be degenerate) in order to make apparent that we are considering the flavor-nonsinglet case. Hence the twist-2 operators do not mix and are multiplicatively renormalizable.
Working with Wilson fermions it is straightforward to write down lattice versions of the above operators. One simply replaces the continuum covariant derivative by its lattice analogue. However, O(4) being restricted to its finite subgroup H(4) (the hypercubic group) on the lattice, the constraints imposed by space-time symmetry are less stringent than in the continuum and the possibilities for mixing increase capi (); roma (); grouptheory (); pertz ().
While the H(4) classification for operators and with has been treated in detail in Ref. grouptheory (), we have to refer to Ref. tensor () for the operators . Note however that the classification of the latter operators for can be derived from the results presented in Ref. grouptheory ().
In our investigations of hadronic matrix elements we have considered the following operators whose renormalization factors have already been studied in Ref. reno ():
Their labels refer to the structure function moments that they determine. These operators have been selected such that they have a definite transformation behavior under H(4) (i.e., belong to an irreducible multiplet) with as little mixing as possible. Moreover we have tried to minimize the number of nonzero momentum components required in the evaluation of the hadronic matrix elements. Note, however, that in the numerical simulations reported in Ref. timid () different operators (though from the same H(4) multiplets) have been used for the matrix elements and .
The operators and transform according to inequivalent representations of H(4), although they belong to the same irreducible O(4) multiplet in the continuum. Therefore their renormalization factors calculated on the lattice need not coincide. The same remark applies to and .
Since our matrix element calculations now involve additional operators, not considered in Ref. reno (), we have extended the above list by the following operators, again guided by the H(4) classification given in Refs. grouptheory (); tensor ():
The operator yields the same structure function moment as . However, in contrast to it cannot mix with any operator of the same or lower dimension. On the other hand, it has the disadvantage that one needs spatial momenta with two nonvanishing components in order to extract the moment from its matrix elements. The latter fact is the reason why we did not employ it in our previous investigations of nucleon structure. The operators constructed from are relevant for transversity.
Furthermore, we have studied the following operators without derivatives (“currents”):
where all quark fields are taken at the same lattice point. Finally we have also considered the quark wave function renormalization constant .
In Table 1 we list all operators studied, along with the H(4) representation they belong to and their charge conjugation parity .
While in the evaluation of hadronic matrix elements the members of a given operator multiplet require different momentum components to be nonzero and hence are of different usefulness, such distinctions do not matter in our computation of renormalization factors. Therefore we consider not only individual operators but also complete operator bases for the representations under study. The representations studied and the chosen bases are given in Appendix A.
Concerning the mixing properties a few remarks are in order. Mixing with operators of equal or lower dimension is excluded for the operators , , , , , , , , , , as well as for the currents.
The case of the operator , for which there are two further operators with the same dimension and the same transformation behavior, is discussed in Refs. grouptheory (); pertz (). Similarly, could mix with another operator of the same dimension. The operators , , on the other hand, could in principle mix not only with operators of the same dimension but also with an operator of one dimension less constructed from . A few more details on the mixing issue can be found in Ref. timid (), in particular in Appendix B.
Our analysis ignores mixing completely. This seems to be justified for . Here a perturbative calculation gives a rather small mixing coefficient for one of the mixing operators roma (); pertz (), whereas the other candidate for mixing does not appear at all in a one-loop calculation of quark matrix elements at momentum transfer zero, because its Born term vanishes in forward direction. The same is true for all operators of dimension less or equal to 6 which transform identically to : Their Born terms vanish in forward matrix elements, hence they do not show up in a one-loop calculation at vanishing momentum transfer. In the case of , however, the mixing with an operator of lower dimension is already visible at the one-loop level even in forward direction. Nevertheless, the nucleon matrix elements of the operators mixing with and seem to be small, at least in the quenched approximation timid ().
Iii The method
We calculate our renormalization constants with the help of the procedure proposed by Martinelli et al. marti () (the Rome-Southampton approach). It follows closely the definitions used in (continuum) perturbation theory. We work on a lattice of spacing and volume in Euclidean space. For a fixed gauge let
denote the nonamputated quark-quark Green function with one insertion of the operator at momentum zero. It is to be considered as a matrix in the combined color and Dirac space. The corresponding vertex function (or amputated Green function) is given by
where for or
denotes the quark propagator. We define the renormalized vertex function by
and fix the renormalization constant by imposing the renormalization condition
at , where is the renormalization scale. So can be calculated from the relation
with . Here is the Born term in the vertex function of computed on the lattice, and denotes the quark field renormalization constant. The latter is taken as
again at . Aiming at a mass-independent renormalization scheme we finally have to extrapolate the resulting values of to the chiral limit.
Note that there will be no lattice artefacts in Eq. (29), because they come with operators of opposite chirality in the vertex function, and these drop out when the trace is taken. Still, matrix elements of the renormalized operators will in general have lattice artefacts because the operators are not improved. Once improved operators are available one can evaluate their renormalization factors using the methods described in this paper.
Equations (28) and (30) (in the chiral limit) together define a renormalization scheme of the momentum subtraction type which is called scheme marti (). Here RI stands for “regularization independent”. This nomenclature refers to the fact that the definition of the scheme does not depend on a particular regularization – here we have used a lattice cutoff just for definiteness and because the lattice regularization will be the basis of our numerical investigations. The scheme, on the other hand, can only be defined within dimensional regularization and is therefore restricted to perturbation theory.
The scheme differs from the RI-MOM scheme only in the definition of the quark field renormalization constant, which in the scheme is more suitable for the numerical evaluation.
In general, the scheme will not agree with any of the momentum subtraction schemes used in continuum perturbation theory. It is therefore desirable to convert our results to a more popular scheme like the scheme. Another reason for converting to the scheme lies in the fact that many of the operators discussed in this paper appear in the operator product expansion along with the corresponding Wilson coefficients, which are generally given in the scheme. Hence we have to perform a finite renormalization leading us from the scheme to the scheme if we want to use our renormalized operators together with the perturbative Wilson coefficients. This finite renormalization factor can be computed in continuum perturbation theory using, e.g., dimensional regularization. The details needed for the evaluation of this factor will be discussed in Sec. VII.
If the operator under study belongs to an O(4) multiplet of dimension greater than 1, i.e., if it carries at least one space-time index, the trace in Eq. (29) will in general depend on the direction of . This has the immediate consequence that the renormalization condition (28) violates O(4) covariance even in the continuum limit. In the continuum, one can restore O(4) covariance by a summation over the members of the O(4) multiplet. On the lattice, each operator when renormalized according to Eq. (29) has in general its own factor, and only after conversion to a covariant scheme all operators in an irreducible H(4) multiplet will have the same renormalization factor. However, it is also possible to define a common factor for all members of such an H(4) multiplet already in the RI-MOM framework by taking a suitable average. If labels the members of the chosen basis of the multiplet we can average over this basis and calculate from
with . This procedure has two advantages. It is simpler than working with a different renormalization factor for every single operator, and it leads to a smoother dependence of the results on , because it reduces the amount of O(4) violation.
The bases actually used in our calculations are given in Appendix A. Whenever we want to refer to this averaging procedure we shall write the respective operator with a bar on top, i.e., means precisely the operator (4) while refers to a result for the operator multiplet (93), and analogously for the currents.
Ideally, the scale at which our renormalization constants are defined should satisfy the conditions
on a lattice with linear extent . The inequality should ensure that we can safely use (continuum) perturbation theory to transform our results from one scheme to another. The inequality is supposed to keep discretization effects small. So we have to find a way between the Scylla of difficult to control nonperturbative effects and the Charybdis of lattice artefacts. Whether in a concrete calculation the conditions (32) may be considered as fulfilled remains to be seen.
Let us finally comment on our notation for the renormalization scale. In the case of a general scheme we use the letter , in the case of the scheme we use . We take when dealing with the scheme and in the case of the MOM scheme to be defined below.
Iv Numerical implementation
Let us sketch the main ingredients of our calculational procedure reno (). To simplify the notation we set the lattice spacing in this section. Moreover we suppress Dirac and color indices. In a first step the gauge field configurations are numerically fixed to some convenient gauge, the Landau gauge in our case mandula (). Representing the operator under study in the form
we calculate the nonamputated Green function (24) as the gauge field average of the quantity
constructed from the quark propagator on the same gauge field configuration. Working in the limit of exact isospin invariance we do not have to distinguish between and propagators. With the help of the relation
we rewrite as
appearing in this expression can be calculated by solving the lattice Dirac equation with a momentum source:
Here represents the fermion matrix. So the number of required matrix inversions is proportional to the number of momenta considered. But the quark propagators, which we need for the amputation and the computation of the quark wave function renormalization, are immediately obtained from the quantities already calculated.
Strictly speaking, one should evaluate the quark propagators going into the calculation of in Eq. (26) on configurations that are statistically independent of those used for the computation of the Green functions (24) in order to avoid unwanted correlations 222We thank F. Niedermayer for drawing our attention to this point.. If we calculate these expectation values on two independent ensembles, statistical fluctuations in the quark propagators and the Green functions (24) are uncorrelated, and (31) gives a good estimate of the true . If we calculate both expectation values on the same ensemble, the fluctuations will be correlated, particularly if the configuration number is small. It can be shown efron () that this will give a bias proportional to . In our case statistical fluctuations are very small, because of our use of momentum sources, so we do not expect a large problem. Indeed, estimating the bias introduced by our procedure (see, e.g., Refs. bias (); efron () for appropriate methods) we confirm this expectation. Note that such correlations exist also in the calculation of hadronic matrix elements from ratios of correlation functions. However, given the typical number of configurations used in these studies, the bias should be very small.
Another computational strategy would be to choose a particular location for the operator. Translational invariance ensures that this will give the same expectation value after averaging over all gauge field configurations. For this method we need to solve the Dirac equation with a point source at the location of the operator and (in the case of extended operators) for a small number of point sources in the immediate neighbourhood. For operators with a small number of derivatives the point source method would require fewer inversions, but it turns out that relying on translational invariance increases the statistical errors.
The required gauge fixing necessarily raises the question of the influence of Gribov copies. Fortunately, investigations of this problem indicate that the fluctuations induced by the Gribov copies are not overwhelmingly large and may be less important than the ordinary statistical fluctuations gribov1 (); gribov2 () (see also Ref. zci ()).
Since the numerical effort is proportional to the number of momenta, the proper choice of the momenta considered is of particular importance. In order to minimize cut-off effects we choose them close to the diagonal of the Brillouin zone and achieve for most operators an essentially smooth dependence on the renormalization scale . It goes without saying that in this way we cannot eliminate lattice artefacts completely. However, more sophisticated strategies for coping with the cut-off effects such as those suggested in Refs. roiesnel (); boucaud () would require the use of far more momenta than we can afford when working with momentum sources. As the treatment of lattice artefacts is a subtle issue anyway we have decided to keep the advantage of small statistical errors provided by the above procedure and to deal with the discretization errors in a different manner.
A specific lattice artefact is caused by the chiral symmetry breaking term of the quark propagator. In position space this term is concentrated at very short distances, for fermion actions obeying the Ginsparg-Wilson condition it is even exactly a delta function. In momentum space it gives rise to the Wilson mass term, in the inverse propagator. The authors of Ref. ferenc () discuss one method of suppressing this artefact (for an earlier treatment of the same effect see Ref. offshell ()). Here, our approach to this problem is to define from Eq. (30), in which the trace removes the Wilson mass term, and to suppress the remaining lattice artefacts by using the perturbative subtraction scheme described in Sec. IX.
V Monte Carlo ensembles
We work with two degenerate flavors of nonperturbatively improved Wilson fermions (clover fermions). For the explicit form of the fermionic action see, e.g., Ref. offshell (). As our gauge field action we take Wilson’s plaquette action. In Table 2 we collect the parameters of our simulations, , , the clover coefficient and the lattice volume along with , the pion mass in lattice units. Table 3 contains the critical hopping parameters as well as the values of the Sommer parameter and the average plaquette in the chiral limit lambda () which are employed in this paper. Note that the results for the chiral extrapolation of given here are based on a larger set of data than that used in Ref. lambda ().
The statistical errors will be calculated by means of the jackknife procedure.
Vi Chiral extrapolation
As already mentioned in Sec. III we have to extrapolate our results obtained at nonvanishing quark masses to the chiral limit. This will be done for each at fixed values of . In the cases where the simulations for different values of have been performed on different volumes, i.e., for and , the sets of momenta used depend on , and some kind of interpolation is required. For this purpose we fit the data on the larger lattices () with cubic splines in . Except for very small momenta, which will not influence the final results, these fits yield a very good description of our data. An example is shown in Fig. 1.
Of course, “wiggles” in the data (caused by lattice artefacts) will be smoothed out by this interpolation. These wiggles are less pronounced on the larger lattices than on the smaller ones. That is why we have chosen to work with the momenta coming from the smaller lattices so that we can use the data on these lattices directly without any interpolation and have to interpolate only the results obtained on the larger lattices.
Alternatively, one could use the interpolation for all values. This however leads to negligible differences in the final results.
For the chosen momenta we can then extrapolate our data to the chiral limit. This is done linearly in , i.e., by means of a fit of the form
where the fit parameter is identified with the desired value of the renormalization factor in the chiral limit. Note that this is essentially a linear fit in the quark mass. The ansatz is motivated by the fact that in perturbation theory the leading quark mass dependence is linear as the chiral limit is approached (see, e.g., Ref. offshell ()).
in order to get an idea of the impact of the chiral extrapolation on the final results. (Note that this corresponds to a three-parameter fit at three data points for .)
However, there is an exceptional case where these simple extrapolations are not trustworthy. This is the pseudoscalar density . In this case one expects that vanishes with the quark mass , because develops a pole in . Therefore we follow Ref. cudell () and try to subtract the pole contribution using a fit of the form
Here is the critical hopping parameter defined for fixed by the vanishing of the pseudoscalar mass . The values used in this paper can be found in Table 3. The fit parameter is then identified with the inverse of in the chiral limit. Examples of such fits are shown in Fig. 4. The curvature in the data is clearly visible establishing the existence of the Goldstone pole.
How can we judge the reliability of the resulting numbers? From the operator product expansion lane (); pagels () we expect that is inversely proportional to , i.e., should become independent of . Therefore we plot versus in Fig. 5. The independence seems to be satisfied with reasonable accuracy. Thus we are confident that our extrapolation for works fairly well. Nevertheless, the results for must be considered with some caution.
Vii Input from continuum perturbation theory
In Sec. III we have explained how one can compute nonperturbative renormalization factors leading us from the bare lattice operators (at lattice spacing ) to renormalized operators in the scheme (renormalization scale ). In this section we collect results from continuum perturbation theory which will be needed for the conversion to standard renormalization schemes such as the scheme.
If the operator under study is multiplicatively renormalizable the operator renormalized in some scheme at the scale is related to the bare lattice operator by
The scale dependence of the renormalized operator is determined by the anomalous dimension
Here the derivative is to be taken at fixed bare parameters, and it is implicitly assumed that the cutoff has been removed in the end. In perturbation theory is expanded in powers of some renormalized coupling constant :
Note that the one-loop coefficient is scheme independent.
Similarly we define the quark field renormalization constant in the scheme so that the renormalized quark propagator is given by . In the scheme, is then specified by Eq. (30) or its continuum analogue. For the anomalous dimension of the quark field we adopt the definition
The running of the coupling constant as the scale is varied is controlled by the function
Again, the derivative is to be taken at fixed bare parameters and it is implicitly assumed that the cutoff has been removed in the end. The perturbative expansion of the function can be written as
In this case the first two coefficients and are scheme independent.
Integrating Eq. (47) we obtain
with the parameter appearing as an integration constant.
In the same spirit we define the so-called RGI (renormalization group invariant) operator, which is independent of scale and scheme, by
where depends only on (or on the bare coupling parameter ). Once we know (or equivalently ), multiplication with will allow us to evaluate and hence the operator (or rather its matrix elements) in any scheme and at any scale we like, provided we know the and functions sufficiently well.
In the two-loop approximation, i.e., setting for , one can evaluate the integral in Eq. (51) easily:
With the help of the methods described in Secs. III and IV we can compute numerically for some range of scales . Knowledge of would then permit us to compute . However, being not covariant for most operators, the scheme is not very suitable for evaluating anomalous dimensions. Therefore we will adopt a two-step procedure for computing . In the first step we transform the numerical results for to a covariant “intermediate” scheme , e.g., the scheme, and in the second step we use the anomalous dimension and the function in this scheme to compute and hence . Thus we could in principle compute as
where denotes the finite renormalization factor leading from the scheme to the scheme and all the scales have been identified with the scale initially set in the scheme.
The most obvious choice for is of course the scheme. However, it will turn out to be advantageous to consider also a kind of (perturbative) momentum subtraction scheme, which we call MOM scheme. This is defined by requiring that in the renormalized vertex function the coefficient of the tree-level (or Born) term equals one at the renormalization scale . To make this definition unambiguous one has to specify in each case the basis used for the other contributions. The quark field renormalization constant is taken to be the same as in the scheme:
The MOM scheme is covariant and rather “close” to the scheme so that the conversion factor from to MOM usually differs less from one than the factor leading from to .
We can expand the conversion factor leading from to in powers of the coupling constant:
we obtain for the conversion factor from to MOM:
If the vertex function is proportional to the Born term (as it happens in simple cases), the scheme and the MOM scheme do not differ. So we have
and consequently . However, in the generic case the matrix will contain also contributions that are not a multiple of . As we consider only operators which are multiplicatively renormalizable in the continuum these additional contributions are finite, but they make different from and are responsible for the dependence of on the direction of the momentum . An explicit example is discussed in Appendix B.
Working with the MOM scheme it is quite natural to expand in a coupling constant which is similarly defined through a momentum subtraction procedure. Therefore we have also considered the scheme as defined in Ref. cheret2 (). The corresponding coupling constant is related to the coupling constant by
where in the Landau gauge