Non-perturbative renormalisation and improvement of the local vector current for quenched and unquenched Wilson fermions

By considering the local vector current between nucleon states and imposing charge conservation, we determine its renormalisation constant and quark mass improvement coefficient for Symanzik $O(a)$ improved Wilson fermions. The computation is first performed for quenched fermions (and for completeness also with unimproved fermions) and compared against known results. The two-flavour unquenched case is then considered.


Introduction
A naive discretisation of fermions onto a hypercubical lattice gives for the action errors of O(a 2 ). However returning to the continuum then leads to the famous fermion doubling problem when we find 15 extra copies of our original fermion. To cure this Wilson [1] added the 'Wilson term' to the action so that the copies decouple in the continuum limit: but then discretisation errors are O(a). As the gluon part of the action (sum of plaquettes in this article) has already O(a 2 ) errors, it is desirable to also achieve this for the fermion action. The Symanzik programme 1 allows a systematic reduction of errors 2 to O(a 2 ). An additional operator of dimension 5 (the 'clover term') is added to the Lagrangian with a coefficient c sw suitably adjusted so that on-shell quantities such as masses now have O(a 2 ) errors. For matrix elements it is also necessary to add further higher dimensional operators to the original operator to achieve O(a) improvement. In this letter we shall be concerned with the determination of the improvement coefficients of the associated improvement operators for the local vector current: V (q) µ ≡qγ µ q. For this operator just two additional operators am q V (q) µ and ia∂ ν T (q) µν are required giving the O(a) improved vector current V (q)IMP µ , and renormalised current V (q)R µ , as with T (q) µν =qσ µν q, σ µν = i[γ µ , γ ν ]/2, and ∂ µ φ(x) ≡ [φ(x +μ) − φ(x −μ)]/(2a). The ia∂ ν T (q) µν operator only plays a role in non-forward matrix elements and will not be considered further here. Thus to O(a) improve the local vector current we need to determine the improvement coefficient b V (g 0 ) (where g 0 is the bare coupling constant). Also as this current is not conserved on the lattice then, as discussed in the next section, it is renormalised with renormalisation constant Z V (g 0 ). Perturbatively we have [7,8,9,10] to one loop (independently of the presence of fermions), where for unimproved fermions c sw = 0 and for O(a) improved fermions However in regions where numerical simulations are performed, β ≡ 6/g 2 0 ∼ 6.0 − 6.4 for quenched fermions (n f = 0) and 5.2 − 5.3 for unquenched fermions (n f = 2), the above formulae may not be applicable.
In this article we shall determine Z V and b V non-perturbatively by considering (nucleon) matrix elements of the time component of the local vector current and imposing charge conservation (the details will be given in the following section). We shall first consider quenched fermions and compare our results with known results in the literature. (While most of the results will be for O(a) improved fermions, for comparison, we shall also briefly consider the unimproved case, c sw = 0.) This is then followed by the unquenched case. Preliminary results have appeared in [11].

The Conserved and Local Vector Currents
There is an exact global symmetry of the lattice action q → e −iαq q,q → e iαqq . This global symmetry is flavour conservation (if you just rotate the quarks of one flavour) or baryon number conservation (if you rotate all quark flavours equally). Separate quark transformations are possible because in pure QCD there are no flavour changing currents. Upon using the Noether theorem this symmetry gives an exactly conserved vector current or CVC of (Being conserved this current requires no renormalisation constant and is O(a) improved 3 .) By this we mean that the Ward Identity, WI, is where Ω is an arbitrary functional of the U, q andq fields. ∆ µ is the backward (Although in eq. (5) we integrate out the fermion fields, . . . q , and take the average over the gauge fields, . . . U , the equation is already true configuration-by-configuration. ) We immediately see that if the region over which Ω is defined does not contain x the RHS of eq. (5) vanishes so that ∆ µ J (q) µ (x + 1 2μ ) = 0. If Ω contains x then the RHS effectively 'counts' the number of quarks and anti-quarks in Ω.
Here, in this study, we shall take Ω → B(t)B(0) where B is the standard (stationary, or p = 0) nucleon operator, containing two u quarks and one d quark summed over the spatial planes 4 . (This was previously used as part of a project to determine moments of structure functions; other choices are of course possible.) µν (x +μ) . 4 We actually take Γ unpol αβ B β (t)B α (0) where Γ unpol ≡ 1 2 (1 + γ 4 ) projects out the unpolarised component of the nucleon field, Inserting this Ω in eq. (5) and also summing over spatial points of where R is defined as the ratio of three-point to two-point correlation functions, and for an N 3 S × N T lattice with c (q) i (i = 1, 2) being constant eq. (6) can be solved to give the result c So R should be a constant with jump, or discontinuity, given by χ (q) . Note that this should be true to 'machine accuracy'. (Indeed in this special case R may also be taken to be the ratio of three-to two-point correlators for a single configuration.) This result may also be shown using transfer matrix methods, as indicated in the appendix, where it is also demonstrated that c 1 . The lattice computation of the three-and two-point functions for R[O(τ )] follows the standard way, see eg [12]. (The sourceB(0) and sink B(t) have been additionally improved by non-relativistic projection and Jacobi smearing to increase the overlap with the ground state nucleon, as described for example in [13]. This does not affect the arguments given above.) We must, in principle, compute a quark-line connected contribution (to the operator) and a quark-line disconnected term as shown in Fig. 1. This latter term is numerically extremely difficult to compute, due to ultra-violet fluctuations. However for the CVC this term is in fact constant. This may be easily seen by substituting Ω → B(t)B(0) q into the WI, eq. (5). The RHS is then zero and so using the same argument as before the appropriate ratio is constant for all τ . There is thus no contribution to the discontinuity. Indeed on a finite lattice, we would expect this constant to be exponentially small (with exponent proportional to N T ). Physically there is no quark-line disconnected term because creating a quark-antiquark pair cannot change the charge. Nevertheless with an eye on the computation of the local current it is useful to consider the difference between the u and d operators, ie the non-singlet operator, in which the quark-line-disconnected terms cancel and so we find where the discontinuity between the two constants, ∆R[J ], in eq. (11) is given by An example for the ratio R[J a very good signal is observed. For the jump, as expected, we find 1 to within machine precision.
The local vector current (LVC) does not obey the WI given in eq. (5), but as there is an additional term formally of O(a) on the RHS of this equation. However perturbatively expanding eq. (5) gives loop graphs with ultra-violet divergences ∼ 1/a and so this additional term gives a contribution. Thus to obey the WI to O(a 2 ), V µ must be renormalised with renormalisation constant Z V and the quark-line disconnected term in Fig. 1 may give an O(a) contribution. Again to avoid computing this term we consider the non-singlet operator 5 . The necessity for a renormalisation constant is also reflected in the fact that ∆R[V ] will not be equal to one. This is illustrated in the RH picture of Fig. 2. So we can define the renormalisation and improvement constants (Z V and b V respectively) by demanding that the local current has the same behaviour as the conserved current, ie 6 5 Tests from computing R[V (q) 4 (τ )] for q = u, d separately showed indirectly that the contribution from the quark-line disconnected term must numerically be very small. So the difference between the singlet and non-singlet renormalisation constant must also be very small. 6 It is more precise to define ]) −1 against am q , we see that the intercept gives Z V while the gradient gives Z V b V . Note that Z V has potential O(a 2 ) differences to other definitions of the renormalisation constant while b V has possible O(a) differences.
While this procedure is correct in the quenched case, for the unquenched case for complete O(a) cancellation, g 0 in Z V should be replaced byg 0 where we haveg 2 0 = g 2 0 (1 + b g am q ), [14]. Thus the LHS of eq. (13) should now become ). So while the intercept still gives Z V , the gradient and hence b V is modified. b g is only known to one loop perturbation theory, b g = 0.01200n f g 2 0 + O(g 4 0 ), [14]. Estimating ∂Z V /∂g 0 for unquenched fermions from the Padé fit, eq. (17) and Table 4 to be ∼ −1.3 gives roughly for the additonal term a decrease of ∼ 0.02. As numerically b V will turn out to be ∼ 2, then this could give a 1 − 2% correction. At present, due to uncertainties in estimating this extra term, we shall ignore this small correction factor.   ]) −1 . The fit intervals chosen were [n l , n u ] and [N T − n l , N T − n u ] where [n l , n u ] = [9,13], [6,11], [7,16] for β = 6.0, 6.2 and 6.4 respectively.

Simulation parameters and raw results for O(a) improved fermions
We have made runs at β = 6.0, 6. ]) −1 are given in Tables 1 and 2. The bare quark mass is defined as where κ is the hopping parameter of the simulation. It is necessary to find κ c (g 0 ), ie the critical point where the (bare) quark mass vanishes. From PCAC we know that the quark mass is ∝ m 2 ps and we can use this to determine where the quark mass vanishes. However a more precise/stable determination was often possible if the PCAC quark mass m AW I q was used, so we fitted the dimensionless quantity, where X = r 0 m AW I q and m IM P q = m q (1+b m am q ) is the O(a) improved quark mass. However determining am AW I q requires a knowledge of the c A improvement coefficent. While this is known for quenched fermions [15]   ]) −1 . The fit intervals chosen were [9,13] and (9) 0.7799 (7) Table 3: κ c (g 0 ), the intercept Z V , gradient Z V b V and b V for the quenched data sets (β = 6.0, 6.2 and 6.4) and for the unquenched data sets (β = 5.20, 5.25 and 5.29). this more sensitively affects the value of am AW I q . Thus in this case we have set X = (r 0 m ps ) 2 . r 0 in eq. (15) is the 'force' scale, [16]. While r 0 could be omitted for quenched fermions, for unquenched fermions r 0 becomes quark mass as well as coupling constant dependent and we use the results given in [17], supplemented by [18]. For quenched fermions, we fitted both B 1 and B 2 coefficients, which practically meant that no b m coefficient was needed, being absorbed into the B 2 coefficient, while for unquenched fermions, as we had only three quark mass values we set B 2 = 0 and used a tadpole improved value of b m , following the prescription given in [10] (and using the plaquette values given in [17]). This gave the results for κ c in Table 3. For quenched fermions these numbers are in good agreement with those given in [15].
As discussed in section 2 we now measure the intercept and gradient of (∆R[V ]) −1 against am q . The results are shown in Fig. 3 for both the quenched and unquenched data. Very good linearity is seen in the results, which enables a precise estimation of the intercept and slope. The results are also given in Table 3  may be more problematic. To give an estimate of possible further errors, if κ c is varied by δκ c , then we have For example, taking δκ c ∼ 0.2 × 10 −4 gives δ(am q ) ∼ 0.0005 and δZ V ∼ 0.0006, δb V ∼ 0.001. These are of the same order or less than the given statistical error.

Z V and b V for O(a) improved fermions
We are now in a position to give our final results for Z V and b V . From eq. (13), dividing the gradient in Fig. 3 by the intercept gives in Table 3 in the third and fifth columns Z V and b V respectively. As for both quenched and unquenched fermions we have three β values, we Unquenched (n f = 2) Z V -0.796 0.0652 -0.667 b V -0.614 -0.0432 -0.767 can attempt to make a Padé-type fit of the form for Z V (g 0 ), b V (g 0 ) constrained to reproduce the weak coupling results, eqs. (2), (3). There are thus 2 free parameters. In Table 4 we give the results of the fits. Note however that even though we have matched to the weak coupling results, as the range of β where the numerical results lie is rather small and far away from this region, intermediate regions may not be represented so well. We estimate total errors on these Padé results to be about 1 2 % for Z V for both the quenched and unquenched cases and for b V , 1% and 2% for quenched and unquenched fermions respectively. (Note also that for the unquenched results there is a further error in b V of 1 − 2% due to an uncertainty in b g as discussed previously in section 2.) Alternative non-perturbative determinations for quenched O(a) improved fermions have been given by the ALPHA Collaboration, using the Schrödinger functional method, [19], the LANL Collaboration, [20,21] using other Ward identities and the SPQcdR Collaboration [22] using the RI-MOM scheme, [23]. All these other methods have reasonable agreement with the results given here, as can be seen in the following pictures. (One should remember that Z V definitions in particular RI-MOM can vary by at least O(a 2 ) and b V definitions can vary by O(a).) In Figs. 4, 5 we plot our results for Z V and b V respectively for both quenched and unquenched fermions. As, in particular for unquenched fermions the available β-range is rather narrow and far from the continuum limit the Padé extrapolation should be treated with some caution. However it is encouraging to note that the ordering of the points is correct. Thus the Padé results should be regarded mainly as an interpolating formula between the various β values. Note however in the quenched case that the Padé results track quite well the ALPHA results.
Some of our earlier results for Z V for quenched fermions can be found in [24]. Note also an example of the large increase in the statistical fluctuations shown there when matrix elements with moving nucleons ( p = 0) are considered.  [19], LANL, [21] and SPQcdR, [22]. Padé fits are also given for our (full lines) and the ALPHA (dot-dashed line) results. The dashed line is the first order perturbative result, eq. (2).

Z V for unimproved fermions
Although most of our results are for O(a) improved fermions we also have results at β = 6.0 for quenched unimproved fermions (c sw = 0) which serve as a useful check on the improved results. The method is identical to that described previously, so we shall just give the results here. In Table 5 we give the parameters of the runs and raw results. Note that the range of quark masses is greater than in the O(a) improved case and that we have runs on different volumes at the same quark mass. It is apparent from the table that finite volume effects are rather small. Plotting ]) −1 against am q gives the result Z V = 0.6436(3) at β = 6.0 (for completeness the gradient is 0.9088 (27)), using κ c = 0.157129 (10). We see that the numerical value is somewhat lower than the equivalent O(a) improved value. This result is consistent, but a little larger than the recent determination given in [25] (although there different β-values are used). Again this is not unexpected, as different determinations of Z V can now differ by O(a) terms (see also [24]).

Conclusions
Our method is not in disagreement with results of other approaches for O(a) improved quenched fermions. For example even for b V there is only a few percent scatter in comparison with other methods. This is what one would expect with O(a) discrepancy between the various definitions.
Z V and b V are both further away from 1 in the unquenched case than in the quenched case at comparable lattice spacings (roughly a n f =2 (5.25) ∼ a n f =0 (6.0)). This is partially an effect from the use of (numerically) smaller β factors and also partially because of the presence of fermions which induces a shift in the effective β. Finally we note that we cannot use one loop perturbation expansion results in the present day quenched/unquenched β regions, although this can be improved by using TI-RGI-BPT, eg [26], or BLM/TI methods, [27].  (where the second equation assumes that the excited states have died away and takes the charge on the nucleon, χ N * being the parity partner of the nucleon and where χ (q) N . All of these expressions are independent of τ . The transfer matrix approach thus again predicts perfect plateaus for J (q) 4 . Furthermore, the discontinuity is given by the difference between these two expressions. As i|B|j is only non-zero if the difference in the charges of |i and |j is equal to that of the charge of a nucleon, then we can replace χ From eqs. (19), (20) we see that for t < N T /2 (the case considered here) the exponential factor is larger for the term 0 < τ < t which means that c