Abstract
We propose a nonperturbative method for computing the renormalization constants of generic composite operators. This method is intended to reduce some systematic errors, which are present when one tries to obtain physical predictions from the matrix elements of lattice operators. We also present the results of a calculation of the renormalization constants of several twofermion operators, obtained, with our method, by numerical simulation of , on a lattice, at . The results of this simulation are encouraging, and further applications to fourfermion operators and to the heavy quark effective theory are proposed.
CERNTH.7342/94
LPTHE prep. 94/52
ROME prep. 94/1022
SHEP 94/9503
A GENERAL METHOD FOR NONPERTURBATIVE
RENORMALIZATION OF LATTICE OPERATORS
G. Martinelli, C. Pittori, C.T. Sachrajda, M. Testa and A. Vladikas
Dip. di Fisica, Università degli Studi di Roma “La Sapienza” and
INFN, Sezione di Roma, P.le A. Moro 2, 00185 Rome, Italy.
Theory Division, CERN, 1211 Geneva 23, Switzerland.
L.P.T.H.E., Université de Paris Sud, Centre d’Orsay,
91405 Orsay, France.
Dep. of Physics, University of Southampton,
Southampton SO17 1BJ, U.K.
Dip. di Fisica, Università di Roma “ Tor Vergata” and
INFN, Sezione di Roma II,
Via della Ricerca Scientifica 1, 00133 Rome, Italy.
CERNTH.7342/94
November 1994
1 Introduction
Renormalization of lattice operators is a necessary ingredient to obtain many physical results from numerical simulations, such as meson decay constants, form factors, structure functions, mixing amplitudes, etc. In this paper we study the renormalization of composite operators in lattice . We propose a method which avoids completely the use of lattice perturbation theory and allows a nonperturbative determination of the renormalization constants of any composite operator. For illustrative purposes, we consider some specific applications to matrix elements of twofermion operators. The approach is particularly useful in those cases, such as the scalar or pseudoscalar densities, where it is not possible to use the chiral Ward identities to determine the renormalization constants nonperturbatively^{1}^{1}1In the following we will denote the method to determine the renormalization constants, using the Ward identities, as the Ward identity method. [1]–[3]. Moreover, our proposal can be applied to many other cases, such as the renormalization of the fourfermion operators of the effective weak Hamiltonian, and to heavylight currents in the Heavy Quark Effective Theory. Preliminary results have been presented in ref. [4]. A similar attempt, limited to the divergent part of the penguinoperator, can be found in ref. [5].
Renormalized lattice operators must correspond, in the limit of infinite lattice cutoff, to finite chirally covariant operators, which obey the same renormalization conditions as those in the continuum. With just a few exceptions, lattice perturbation theory is used to evaluate the renormalization constants of lattice operators. The problem of mixing with lower dimensional operators requires however, a nonperturbative subtraction of all power divergences [1],[6]–[9]. Apart from this special case, the use of perturbation theory, in the computation of multiplicative constants and dimensionless mixing coefficients, is well justified, provided that the lattice spacing is sufficiently small, i.e. . However, in those cases where the results from perturbation theory, obtained using the bare lattice coupling constant as an expansion parameter, can be checked using nonperturbative methods, there is generally a significant discrepancy. Several solutions to this problem have been proposed so far:

In order to improve the convergence of the perturbative series, the authors of ref. [10] have proposed a “tadpole improved” perturbation theory. Nevertheless, ignorance of higher order contributions still represents an important source of uncertainty in the extraction of physical results.

One may fix nonperturbative renormalization conditions directly on hadronic matrix elements. This procedure was used in ref. [7], to subtract the divergent part of the operator. The price is that one has to sacrifice a physical prediction, for any subtraction imposed on hadronic states. Thus, when several renormalization conditions are necessary, one looses much predictive power [1, 6].
Our proposal is to impose renormalization conditions nonperturbatively, directly on quark and gluon Green functions, in a fixed gauge, with given offshell external states, with large virtualities.
The method consists in mimicking what is usually done in perturbation theory. One fixes the renormalization conditions of a certain operator by imposing that suitable Green functions, computed between external offshell quark and gluon states, in a fixed gauge, coincide with their tree level value. For example, if we consider the generic twoquark operator , we may impose the condition
(1) 
where is one of the Dirac matrices and is the tree level matrix element. The extension to more complicated cases, including fourfermion operators and operator mixing, is straightforward. This procedure defines the same renormalized operators, i.e. the same Wilson coefficient functions, in all regularization schemes, provided they are expressed in terms of the same renormalized coupling constant. However, the coefficient functions now depend on the external states and on the gauge, which must be specified.
In principle this scheme completely solves the problem of large corrections in lattice perturbation theory, which are automatically included in the renormalization constants. Matching of lattice operators to the corresponding ones, renormalized in a continuum renormalization scheme (for definiteness, the scheme), then requires continuum perturbation theory only. When nexttoleading corrections are known, the error on the physical matrix elements is of , in the continuum perturbative expansion, which is expected to have a better convergence. A continuum perturbative calculation of the matching conditions is a common step in all approaches, standard or ‘‘tadpole improved” perturbation theory and nonperturbative renormalization. Our proposal can be applied to any composite operator, unlike the Ward identity method, which is limited to a few, albeit very important, cases^{2}^{2}2It is still true that the statistical accuracy is in general slightly better with the Ward identities method, in those cases where it can be applied, see sec. 6..
The feasibility of using Green functions computed between quark states, in a fixed gauge, was already discussed in ref. [3], in the context of nonperturbative renormalization, using Ward identities. Indeed, in the case of the axial current, it was shown that the signal was rather good and the result was in agreement with the corresponding hadronic gauge invariant determination.
As far as the renormalization conditions are concerned, our proposal is expected to work, whenever it is possible to fix the virtuality of the external states and to satisfy the condition , in order to keep under control both nonperturbative and discretization effects. We stress that this requirement is common to all methods. The results of our numerical investigation are encouraging: they suggest that there does exist a “window” in , where the method can be applied, at values of the lattice coupling constant used in current lattice simulations.
The plan for the remainder of the paper is as follows. In the next section the basic formulae necessary to define our nonperturbative method are presented, using as an example logarithmically divergent twoquark operators, such as the scalar and pseudoscalar densities. The discussion is generalised to all composite operators in section 3, and in the following section the implementation of the method in numerical simulations for twoquark operators is discussed. In section 5 we give some details about the perturbative evaluation of the renormalisation constants on lattices of finite volume. The values of the renormalisation constants obtained in a numerical simulation are presented in section 6, and are compared to the results from perturbation theory. The final section contains our conclusions, a summary of the numerical results and a discussion of possible further applications of the ideas presented in this paper.
2 Nonperturbative renormalization conditions
In this section, the nonperturbative renormalization scheme is presented in detail. We also discuss the conditions which must be satisfied for our method to be applicable.
To facilitate the discussion, it is convenient to classify composite operators into three main classes, according to their ultraviolet behaviour, as :

Finite operators: these include the vector and axial vector currents. Another example is the ratio of the scalar and pseudoscalar renormalization constants, . For these cases, the Ward identity method is also applicable.

Logarithmically divergent operators: this is a large family of operators, relevant to hadron phenomenology. It includes scalar and pseudoscalar densities, fourfermion operators, some components of the energymomentum tensor, etc.

Power divergent operators: these operators are present when mixing with lower dimensional operators is possible. This happens in regularizations which have an intrinsic mass scale, such as the lattice regularization. Examples are often encountered in phenomenological applications of lattice , such as fourfermion operators relevant to transitions or operators of order in the heavy quark effective theory (HQET) [11]–[15].
For simplicity, we first consider the twofermion operators of class 2., and extend the discussion to other operators in this class, and also to those in classes 1. and 3., in the next section. Thus, the formulae given below apply directly to the pseudoscalar and scalar densities. Throughout our discussion, we assume that discretization errors are negligible: in field theory language, this is equivalent to the statement that the renormalized Green functions do not differ appreciably from their values in the limit of an infinite ultraviolet cutoff. The discussion is presented in the limit of small quark masses and in Euclidean spacetime. The discretisation of the quark action is assumed to be performed à la Wilson, characterized by explicit chiral symmetry breaking of the Lagrangian^{3}^{3}3This includes the SWClover action [16].. The extension to staggered fermions is straightforward.
Let us consider the forward amputated Green function , of a twofermion bare lattice operator , computed between offshell quark states with fourmomentum , with , and in a fixed gauge, for example the Landau gauge. Without gaugefixing all Green functions computed between quark and gluon external states are zero. This is also true when the gauge is not completely fixed, e.g. in the Coulomb gauge. We define the renormalized operator , by introducing the renormalization constant
(2) 
is found, by imposing the renormalization condition
(3) 
where is the field renormalization constant, to be defined below. This procedure defines a renormalized operator which is independent of the regularization scheme [17]–[19]. It depends, however, on the external states and on the gauge. This does not affect the final results, which, combined with the continuum calculation of the renormalization conditions, at any given order of perturbation theory, will be gauge invariant and independent of the external states. Let us specify the different quantities entering eq. (3). is defined in terms of the expectation value^{4}^{4}4 “Expectation value” means, as usual, that one averages the Green functions over the gauge field configurations, generated by Monte Carlo simulation, see sec. 4. of the nonamputated Green function , and of the quark propagator
(4) 
where
(5) 
is a suitable projector on the treelevel operator: () for the scalar (pseudoscalar) density. The factor ensures the correct overall normalization of the trace (colourspin=12). Projectors are very convenient when defining Green functions, particularly in the nonperturbative case. They have been extensively used in refs. [18, 19]. Of course one can also use other definitions of .
is the renormalization constant of the fermion field. It can be defined in different ways, some of which are equivalent perturbatively. Beyond perturbation theory, the most natural definition of is obtained from the amputated Green function of the conserved vector current . Indeed, one knows that for the renormalization constant is equal to one:
(6) 
which implies
(7) 
Equations (3)–(7) completely define our method. In the remainder of this section, we discuss some important aspects concerning its applicability.
In the above formulae, for simplicity, we have always considered only forward matrix elements. In general, one has the freedom to define the renormalization conditions at different external momenta and (). The virtualities of the quark states must be much larger than . The reason is that, in order to obtain the physical result, we have to combine the matrix element of the renormalized operator , with a Wilson coefficient function. The latter is computed in continuum perturbation theory, by expanding in at a scale of order . Thus, for the validity of this perturbative calculation, must be large. An important question in our program is whether it is possible to find, on the lattice, a scale which is sufficiently low, in order to have small effects, and sufficiently large, in order to have small higher order corrections. The range of , where the above conditions are satisfied, depends on the value of of the numerical simulation. We stress that, if such a window does not exist, in current lattice simulations, an accurate matching of lattice operators to continuum ones becomes impossible, not only in our approach, but also with any other method.
It might be expected that the condition ensures that perturbation theory is valid. When spontaneous symmetry breaking occurs however, as is the case in , a large value of may not be enough, because of the presence of the Goldstone boson, the pion in our case. At low momentum transfers , Green functions can receive a nonperturbative contribution from the pion pole, which cannot be computed. The naïve expectation is that this effect vanishes as , at large . However, with fermion fields, this contribution is proportional to , even when . The proof is given in the appendix. We conclude that a large value of solves the problem. At low values of , however, for the pseudoscalar density and the axial current, we expect to find visible effects, see secs. 3 and 6.
A danger to the nonperturbative scheme may come from the presence of Gribov copies. We do not address this problem here; it has been investigated in refs. [20, 21]. The results of ref. [21] indicate that, at least for the quantities of interest in the present study, the “Gribov uncertainty”, in current lattice simulations, is at most of the same order as the statistical error.
One may imagine applying the renormalization conditions at very large values of , on lattices which have small physical volumes. This might seem to overcome the problem of the existence of the window in . For physical applications however, it is necessary to evolve the renormalization constants to smaller values of . Nonperturbative contributions, which are not detected at large ’s, may then become important. This demonstrates the necessity, for any method, of the existence of a window , where is the lattice spacing used in numerical simulations of physical quantities. Only then are nonperturbative effects and lattice artefacts small simultaneously.
3 Applications of the method to generic composite operators
In the previous section, using the renormalization of the scalar and pseudoscalar densities as a prototype, the general strategy of the nonperturbative scheme has been presented. In the present section the necessary modifications needed in other cases are briefly described. We will follow the classification of operators introduced in sec. 2.
3.1 Finite operators
In the case of vector and axial vector currents, we do not have the freedom of fixing the renormalization conditions in an arbitrary way: the renormalization conditions are acceptable only when they are compatible with the Ward identities. To be specific, one can apply eq. (3) to the local vector current in order to determine
(8) 
We now show that this definition of , which uses the general prescription given in the previous section, coincides with the one derived from the vector current Ward identity on quark states [1, 22]^{5}^{5}5For the remainder of this subsection we work in lattice units, setting .
(9) 
where is the amputated vertex with momentum transfer . By applying , at , on both sides of eq. (9), one gets
(10) 
The second term on the l.h.s. vanishes when . By tracing eq. (10) with , we then recover the definition of of eq. (8).
One is tempted to follow the same path for the local axial vector current, i.e. to define
(11) 
However this definition of is not equivalent to the one obtained by the Ward identity
(12) 
To see this we proceed as for the vector current, applying to both sides of eq. (12).
(13) 
In the axial case, however, the second term on the l.h.s. does not vanish, as and, indeed, this term is necessary to saturate the Ward identity, in presence of a massless Goldstone boson. As demonstrated in the appendix however, this term becomes negligible for high values of . In this case, by tracing eq. (13), one recovers (11). In conclusion, we expect to find a value of close to those obtained from the Ward identity on hadron or on quark states [3], by using eq. (11), at large .
is another finite quantity, which is completely fixed by the Ward identities, even though the renormalization conditions on and are separately arbitrary. In other words one is free to decide the renormalization conditions of one of the two operators, the scalar density say, and the value of will automatically follow. The procedure of sec. 2, respects the Ward identities and thus guarantees that the correct value of is obtained. This is clearly true only at large values of , for the same reasons as for the axial vector current.
3.2 Logarithmically divergent operators
For this category of operator only the question of operatormixing is left to be discussed. There are two kinds of mixing: the one which can occur also in the continuum and the one induced by the explicit lattice chiral symmetry breaking. We are only concerned here with the latter case, and hence restrict our discussion to operators which renormalize multiplicatively in the continuum. To illustrate our arguments, we consider the fourfermion operator
(14) 
which is relevant to the calculation of the – mixing amplitude. Because of the Wilson term, this operator mixes with operators of different chiralities [1, 6, 23]
(15) 
where the operators are given by , etc. [23]. The mixing coefficients are finite functions of . The overall renormalization is necessary to eliminate the logarithmic divergence.
We now sketch the renormalization procedure for . Let us start from the fourpoint amputated, Green function of , between quark states with momentum
(16) 
where and are spin and colour indices respectively. To find the mixing coefficients , one uses suitable projectors , as was done for twoquark operators^{6}^{6}6 In order to simplify the presentation, we do not discuss here the complications arising from the different possible colour structures., for example
(17) 
By projecting on all the possible colourspin structures, corresponding to the operators , one obtains a system of linear equations in the mixing coefficients . Using the , determined in this way, one can compute the Green function of the subtracted operator . From the amputated vertex of , one can finally remove the logarithmic divergence, by imposing the renormalization condition
(18) 
where
(19) 
and is the suitable projector on the Dirac tensor.
3.3 Power divergent operators
Composite operators are divergent in when they can mix with lower dimensional operators. Examples are given by the operators relevant to deep inelastic scattering and by the penguinoperators appearing in the Hamiltonian. Let us consider the operator
(20) 
mixes with the scalar and pseudoscalar denisties, even at zeroth order in , through the diagram in fig. 1^{7}^{7}7 also mixes with the chromomagnetic operator, e.g. , which is not being considered here.
(21) 
To find the coefficients , one imposes that the twopoint Green function of the operator on the r.h.s. of eq. (21), on offshell quark states, is zero [5]. Once all the divergent parts have been removed, one can impose the renormalization conditions, necessary to remove the overall logarithmic divergence, as discussed above.
4 Implementation of the method in numerical simulations
To implement the ideas of the previous section, let us consider the generic local twofermion operator
(22) 
where is any Dirac matrix. The lattice, bare Green function of , between offshell quark states, can be obtained from
(23) 
where labels the gauge field configurations and is the quark propagator on the single configuration , obtained by inverting the discretized Dirac operator. We denote by the inverse Dirac operator, which is not translationally invariant, in contrast to , which is the usual quark propagator. The gauge field configurations are generated with some standard, gaugeinvariant algorithm; the gauge fields and the quark propagators are then rotated into the Landau gauge, which will be defined in sec. 6. From eq. (23), using , one gets
(24)  
where we have defined
(25) 
because, on a single configuration, the quark propagator is not translationally invariant. Traslational invariance is recovered after averaging over the configurations. We can then write
(26) 
Eqs. (24) and (26) give and in terms of quark propagators, computed by inverting the Dirac operator in the presence of a background gauge field configuration, generated by numerical simulation. With these quantities at hand, one can apply the strategy and the formulae, introduced in sec. 2. We add some information, that can help the reader in practical applications. The quark wavefunction renormalization constant has been computed using
(27) 
where is defined as in eq. (5); the trace is taken over colours and spins and is the renormalization constant of the local vector current . Equations (7) and (27) are two equivalent definitions of . We have used eq. (27) instead of eq. (7), because the former is more convenient when, as in our case, one has only computed quark propagators which stem from the origin. The threepoint () Green function is obtained by inserting the local vector current in the origin. can be obtained with very high accuracy from ratios of matrix elements of the conserved and local currents [2, 24]. From the amputated Green functions
(28) 
by imposing the renormalization conditions (3), at several values of , to be specified in sec. 6, we have determined , as functions of .
If we wish to compute using the procedure defined in sec. 3 we clearly cannot take from eq. (27). From the Ward identity we have
(29) 
where . To avoid derivatives with respect to a discrete variable, we have used
(30) 
In one loop perturbation theory, , in the Landau gauge^{8}^{8}8In general, they differ by a finite term of order .. By using
(31) 
and imposing the renormalization conditions (3), we have determined at several . From the general discussion of the previous section, one expects that , defined through the present procedure, is independent of , up to terms of , as has been effectively found, see sec. 6.
In actual numerical simulations, the renormalization procedure could be obscured by lattice artefacts, i.e. terms of , which are present when one imposes the renormalization conditions at large values of , required to avoid large nonperturbative or higher order effects. As shown in ref.[25] (we use here the same notation), in order to reduce discretization errors, due to the finite value of the lattice spacing, one can compute correlation functions using a nearestneighbour improved fermion action [16]
(32) 
where is the usual Wilson action, is the Wilson parameter and is the lattice field strength tensor. Moreover, it is also necessary to use improved operators defined as
(33) 
where the rotations of the fermion fields are defined as
(34) 
By using the improved action and operators, one eliminates discretization errors of , and only those of or , or higher, remain. In eq. (33), and can be replaced by
(35) 
in computations of the hadron spectrum and of onshell operator matrix elements. As discussed in refs.[26, 27], the improvement procedure is equivalent to expressing all correlation functions in terms of “effective” quark propagators
(36) 
where is the fermion propagator of the improved theory from the point to the origin. In order to eliminate terms of in the offshell correlation functions, used in this study, a further step is necessary [3]: we have to shift the quark propagator by a local term, which in Fourier space appears as a constant, to be added to its diagonal components
In the following, the vector current , the axial vector current , the pseudoscalar density and the scalar density are improved operators, obtained by using eq. (36).
5 Lattice perturbation theory
In order to gain some insight into the nonperturbative results presented in the next section, we have performed the corresponding oneloop perturbative calculation both in the continuum and on the lattice in the Landau gauge. This calculation is an extension of the one presented in ref. [27].
The standard calculation of the oneloop corrections corresponds to a computation of the renormalization constants with a fixed ultraviolet cutoff on the loop momenta, , neglecting terms of in the final result. The renormalization constants obtained in this limit depend on the quark mass and on the external momenta only through the logarithmically divergent terms, which are related to the anomalous dimension of the operator. Since, at one loop, these terms are universal, i.e. they are the same on the lattice and in any continuum regularization, the renormalization constants, necessary to relate lattice to continuum renormalized operators, are independent of the quark mass and external momenta (they depend however on the lattice spacing and on the continuum renormalization scale ). We call this limit “Standard Perturbation Theory” (STP). The nonperturbative calculation however, is performed on a lattice of finite size, and at a fixed value of the lattice spacing and quark mass: thus we cannot make terms of arbitrarily small. In particular, since we want to study the dependence of the renormalization constants, it is important to know at which values of (and ), the perturbative results, obtained on a finite lattice and with a fixed lattice spacing, differ from those obtained in SPT. In order to follow closely the nonperturbative computation, the perturbative results have also been obtained, on a lattice of size . In this way, we keep all terms of or , or higher, in the final result (we call this procedure “Discrete Perturbation Theory”, DPT). The loop integrals were computed numerically, both in SPT and in DPT.
The nonperturbative results have been obtained by tracing suitable amputated Green functions, computed between offshell quark states, in the Landau gauge, cf. eq.(4). In perturbation theory, this procedure differs from the standard one, where one identifies the correction from the term proportional to the treelevel matrix element of the operator under consideration [18, 19]. The difference is due to finite terms which, in the standard procedure, are considered as part of the matrix element of the operator, whereas with the projectors, are considered as part of the oneloop corrections. At oneloop order, these terms are the same in the continuum and on the lattice (in the limit ) and cancel when one computes the difference between the continuum and the lattice renormalization constants.
To be more specific, let us consider the diagrams in fig. 2, computed in Naïve Dimensional Regularization (NDR), in the Feynman gauge
(37) 
where , and is the external momentum. Usually, one takes the term proportional to , i.e. the coefficient of , as the one loop contribution to the bare operator. With the projector (see eq.(4)) the last term in eq. (37) also gives a contribution, so that the factor now becomes . In table 1, we give the oneloop corrections obtained, with the projectors and in the ’t Hooft and Veltman (HV) dimensional renormalization scheme [28], for the twoquark operators considered in this paper. The NDR regularization differs from the HV one by the definition of the anticommuting properties of . For details see for example [19]. In NDR the corrections to the axial current and pseudoscalar density coincide with those to the vector and scalar density respectively. By writing the gluon propagator as , we denote as “Landau” the contribution to the oneloop corrections proportional to .
Operator  Feynman  Landau 

With the procedure, after the subtraction of the pole, the irreducible vertex of fig. 2 has the form^{9}^{9}9 is obtained from using a projector, as described in eq. (4).
(38) 
From table 1, one finds , , etc., in the Feynman gauge and , , etc., for the Landau (longitudinal) terms.
On the lattice, one gets a similar expression
(39) 
The values of , for the Feynman and Landau corrections are given in table 2.
Operator  Feynman  Rotated Feynman  Landau 

In table 2, we have denoted by “Rotated Feynman” the contribution coming from the rotated part of the twofermion operators, see eqs. (33) and (34). The Landau contribution to the rotated part of the operators vanishes. This is required by the fact that the renormalization constants which relate lattice to continuum operators are gauge invariant.
Following refs. [27]–[29], we write
(40) 
where can be computed from tables 1 and 2
(41) 
In the difference between continuum and lattice oneloop corrections, all Landau pieces cancel, as expected. The results for all agree with those of ref. [27].
On a lattice, for a given quark mass and external momenta, there are terms of order which can become important at large values of and . Since with our method we have to impose the renormalization conditions at values of the momenta which are large, it is important to monitor, at least in perturbation theory, the values of masses and momenta at which effects become important. The Green function in DPT has a form similar to (39)
(42) 
where, however, the “constant” now depends on the terms of . Thus, in the calculation of the oneloop corrections in DPT, all the cutoff dependent terms have been included. On a finite volume, the propagators of massless particles at zero momentum must be regulated. The results of Feynman diagrams depend on the regulator, however this dependence vanishes as the volume is increased. In the results presented below, we have set the zeromomentum terms in both the quark and gluon propagators to zero. We observe a small difference in the results from SPT and DPT, at low values of , which is due to this choice. At the value of the quark mass at which we have performed the simulation, mass corrections are negligible.
The coupling constant appearing in eqs. (39) and (42) is the bare lattice coupling. In ref. [10] it was argued that it is possible to optimize the convergence of the perturbative series by reorganizing the expansion in terms of a “boosted” coupling , defined in terms of a physical quantity. In the comparison with the nonperturbative results below we follow the suggestion of ref. [10], and compute the perturbative corrections using an effective coupling
(43) 
where and is the expectation value of the plaquette. We call this expansion “Boosted Discrete Perturbation Theory” (BDPT). Analogously we denote the “Boosted Standard Perturbation Theory” by BSPT.
The results of the calculation of the renormalization constants in BDPT will be compared to the nonperturbative results in the next section.
6 Nonperturbative results
In this section, we give the main results of our numerical study of the nonperturbative renormalization procedure proposed in this paper and a comparison of these results with perturbation theory.
We have performed a simulation, by generating independent gluon field configurations, on a lattice, at . The quark propagators have been computed for a single value of the quark mass (), corresponding to a value of the hopping parameter . All the Green functions have been computed in the lattice Landau gauge, which is obtained by minimizing the functional
(44) 
Possible effects from Gribov copies or spurious solutions [30] have not been considered. As shown below, for those quantities which can be determined also in a gauge invariant way, the nonperturbative results obtained on quark state in a fixed gauge are in good agreement with those obtained using the Ward identity method [3].
Vector current: In fig. 4, the renormalization constant of the local vector current , obtained using eqs. (30) and (31), is given. As expected, is independent of the scale, whithin statistical errors, up to large values of , where distortions due to lattice discretization become important. This is a consequence of the equivalence, which can be established up to terms of , of the method used here to determine , and of the Ward identity for the local vector current, see sec. 3. By using the points at , corresponding to GeV, we get , to be compared with from the Ward identity method, and ( on a lattice, at ) from BSPT. The three methods are in good agreement.
Scalar density.: is an ideal quantity with which to check the validity of our method for a number of reasons: it cannot be determined using the Ward identity method, it is logarithmically divergent and gauge dependent, and it is not affected by the pole of the Goldstone boson (in contrast to the pseudoscalar density and the axial current), sec. 2.
In fig. 5, is shown as a function of the renormalization scale. We report separately the results, obtained by using , eq.(27), and , eq.(30). Notice that the points, corresponding to different definitions of , are in very good agreement where is not too large, and discretization errors are small. This can also be seen by comparing the results in BSPT (dashed curve) to those in BDPT (continous curve). We do not expect agreement with perturbation theory, either at low , where higher order effects become very large, or at high , where lattice distortions are important, as shown by the continous curve. The numerical results follow the theoretical expectations and we find good agreement between the nonperturbative determination of and the predictions from boosted perturbation theory for ( GeV GeV). Notice the () is defined here from the offshell Green function and depends on the anomalous dimension of the scalar (pseudoscalar) density, and hence on log(). Moreover, the result is gauge dependent, even in the continuum, as can be shown by an explicit calculation in oneloop perturbation theory.
Axial vector current: The local axialvector current shares some of the features of the vector current. Its renormalization constant is finite at all orders in perturbation theory [31] and can be determined from Ward identities. However, the axialvector current is coupled to the wouldbe Goldstone boson of , which can give an important contribution at low . , which like , should be independent of , is shown as a function of in fig. 6.
We interpret the strong dependence of , at low , as the non perturbative effect of the pseudoscalar state. Unfortunately, there is no clear sign of the existence of a plateau between the nonperturbative regime and the large region, where lattice artefacts become important^{10}^{10}10The plateau could eventually appear more clearly at larger values of .. Without the information coming from the Ward identity method [3], it would be difficult to determine confidently. Nevertheless, we do observe that, at , just before lattice artefacts become large in DPT, the results are close to the value determined with the Ward identity, [3, 21], see fig. 6. As observed in ref. [3], perturbation theory gives values of smaller than one, while the Ward identity method gives values larger than one, for –. A value larger than one is also suggested by our nonperturbative results.
Pseudoscalar density: The pseudoscalar density shares the main features of the scalar density, with two important differences. Firstly it is coupled to the wouldbe Goldstone boson, and secondly the oneloop perturbative corrections, with the Clover action, are quite large. It is then not surprising that the nonperturbative value of the corresponding renormalization constant lie well below the perturbative result, as shown in fig. 7. Had we used standard perturbation theory in fig. 7, instead of boosted perturbation theory, the discrepancy would have been even worse. The discrepancy could have been anticipated from the results of ref. [3], where it was shown that the ratio , determined by using the Ward identity method, was significantly larger than the result obtained in boosted perturbation theory. Combining this information with the agreement between the nonperturbative determination of and the result from BPT, one would conclude that the difference is due to . This may also be due to the fact that has larger oneloop finite corrections, , and that perhaps the higher order terms, which are not accounted for by one loop boosted perturbation theory, are important.
We also present, in fig. 8, as a function of . The ratio has a behaviour similar to . Some discretization effects appear to be smaller than those observed in and separately. The numerical results support, qualitatively, the theoretical interpretation. In the intermediate range of , this ratio stays almost constant. At small or large values of the results are unstable, due to nonperturbative or discretization effects respectively. At , the value of this ratio obtained from our nonperturbative method is in reasonable agreement with that determined with the Ward identity method, [3]^{11}^{11}11The value published in ref. [3], using half of the present stastistics, was , corresponding to ..
7 Conclusion
We have proposed a new nonperturbative method to renormalize lattice composite operators. It can be used in all cases and is particularly useful when the Ward identity method is not applicable. This method avoids the need to perform any calculations using lattice perturbation theory in the computation of physical quantities from lattice simulations. The success of our proposal is subject to the existence, at current values of , of a window between the nonperturbative region at low momenta and the region of large momenta, where discretization errors become important. This limitation is however, common to all methods (with the exception of the Ward identity method, which can only be applied to a few cases, corresponding to finite operators). The window, necessary to implement the renormalization programme described here, seems to exists already at and we expect that the range of momenta, useful for nonperturbative renormalization, will become larger as increases. We are planning to apply the approach presented in this paper to the renormalization of the four fermion operator (14) and of the heavylight axial vector current in the static theory.
Acknowledgements
We thank S. Petrarca for an early participation to this work and for many discussions. We also thank E. Gabrielli and the members of LPTHE for discussions. G.M. and C.P. thank the theory division at CERN, where part of this study has been performed. We acknowledge the partial support by M.U.R.S.T. and by the EC contract CHRXCT920051. C.P. acknowledges the support by the Human Capital and Mobility Program, contract ERBCHBICT930887. CTS acknowledges the Particle Physics and Astronomy Research Council for its support through the award of a Senior Fellowship.
Appendix
In this appendix we demonstrate that the nonperturbative contributions to forward twoquark matrix elements are suppressed as inverse powers of , where is the external momentum. Let us consider the nonamputated, twofermion Green function
where and is one of the Dirac matrices. It is natural to worry about the infrared behaviour of , because it contains the insertion of a zeromomentum operator. On general grounds, we expect infrared divergences to appear in Green functions computed at exceptional momenta. In physical terms, this is due to the fact that the matrix elements of , at zero momentum, carry information about the physical mass spectrum, which is unaccessible to perturbation theory. In particular, if the operator is the pseudoscalar density () and we work in the chiral limit, a contribution from a Goldstone pole, corresponding to the pion, should appear in . The presence of nonperturbative contributions to is signaled by the appearence of infrared divergences in the perturbative expansion. This is what happens in general. For quark degrees of freedom in the limit of large however, the Operator Product Expansion (OPE) guarantees that can be reliably computed at all orders in perturbation theory. The argument goes as follows. In the limit of large , the dominant contribution to comes from regions of integration (in and ) where the integrand is singular. This implies that, in this limit, we have to study the behaviour of
when and . In both cases of course, we can use the OPE in the following form^{12}^{12}12From a technical point of view, OPE is a weak operator statement and, when is inserted in a Green function containing other operators, a slightly modified form would be appropriate. It is easy to convince oneself, however, that the additionals terms, required in the general case, are not relevant to the present discussion.
for , and
for .
Since QCD is asymptotically free, the behaviour of and can be computed in perturbation theory. The interaction gives rise to logarithmic corrections to the free field behaviour
Using the above equations, in the limit , we find
where
The term is computable in perturbation theory. On dimensional ground, for , it behaves as
As for , we find
is, however, multiplied by , which is a nonperturbative quantity. We have shown however, that, in the large limit, the perturbative contribution dominates by one power of over the nonperturbative one. This is the reason why, even at an exceptional external momentum, is infrared safe in perturbation theory, when is large. It is straightforward to trace the suppression of the infrared divergence, by analyzing the corresponing Feynman diagrams. For example, the infrared divergent part of the diagram in fig. 9a is suppressed by one power of