On cutoff effects in lattice QCD from short to long distances

We discuss kinematical enhancements of cutoff effects at short and intermediate distances. Starting from a pedagogical example with periodic boundary conditions, we switch to the case of the Schroedinger Functional, where the theoretical analysis is checked by precise numerical data with Nf=2 dynamical O(a)-improved Wilson quarks. Finally we present an improved determination of the renormalization of the axial current in that theory.


Introduction
Dynamical fermion simulations with Wilson-type fermions are now possible at small quark masses, large volumes [1][2][3][4][5][6][7][8] and small lattice spacings. With statistical precision reaching the (sub-)percent level, an important uncertainty which remains to be carefully controlled is due to the finite lattice spacing a. In the non-perturbatively O(a)-improved theory, this issue has been investigated in some detail. On the one hand, significant aeffects have been found at lattice spacings of around a = 0.1 fm [9,10]. On the other hand, both in the high precision computations of the scale dependence of coupling, quark masses and other composite operators [11][12][13] (see ref. [14] for a review) and in a recent scaling test [15], such effects were not visible.
In this paper we analyse this apparent contradiction and find that the difference is of a simple kinematical origin, teaching us a more general lesson. In order to keep a-effects small, it is important to ensure that the general conditions necessary for the application of the Symanzik expansion of lattice observables in powers of a are well fulfilled. We discuss the general issue in Sect. 2 and then illustrate it in Sect. 3 on a particular example where precise non-perturbative results are available: a Schrödinger functional correlation function. The lessons we learn there allow us to perform an improved computation of the renormalization of the axial current in Sect. 4, which we check further by considering modified renormalization conditions (Sect. 4.2) before we conclude.
The functions C depend in addition to the kinematical argument x 0 on the intrinsic scale Λ ≡ Λ QCD = O(300 MeV) of the theory as well as the quark masses M i . In the above formula the powers of the lattice spacing, are in principle modified by a further a-dependence in C i (x 0 ) → C i (x 0 , a). It is a consequence of the dependence of the couplings in the Symanzik effective Lagrangian on a as well as the dependence of the coefficients of the effective fields on a. However, since QCD is asymptotically free, these coefficients depend only logarithmically on the lattice spacing when it is small. Such logarithmic terms are of minor importance for our discussion; they are dropped here.

Point-to-point correlators at short distances
Let us assume a large volume and a "point-to-point" correlator, i.e. C lat (x 0 , a) = a 3 x O 1 (0)O 2 (x) with local operators O 1,2 . We focus on the short distance regime, x 0 ≪ 1/Λ , x 0 ≪ 1/M i and an O(a)-improved theory where C 1 ≡ 0. Then x 0 provides the only dimensionful parameter. On purely dimensional ground the leading correction then becomes It is enhanced at small x 0 . In other words what we have used here is that up to logarithmic corrections the short distance correlation functions have a unique powerlaw behaviour in x 0 . A particularly simple example is which one may, in principle, consider for renormalizing the pseudo-scalar density. One expects that the short distance behaviour C PP cont ∼ const. × x −3 0 is difficult to reproduce accurately on a lattice. In fact, one term in the Symanzik expansion originates from an O(a 2 ) correction to the field to eq. (2.1). An order of magnitude enhancement is due to the second derivative of the steep function C cont . Even if c 2 may be arising only at 1-loop of perturbation theory 1 , the considered case suffices to illustrate our main point: cutoff effects may have a significant kinematical enhancement. Particular examples are correlation functions with strong short distance singularities. Of course, this is the reason why the connection between the perturbative short distance regime of QCD and the non-perturbative long distance one is carried out recursively in the strategy of our collaboration [14,19,20]. It is then possible to have a/x 0 ≪ 1 and x 0 Λ ≪ 1. But furthermore, by making use of Schrödinger functional boundary conditions [21], one may easily construct correlators with a weak time-dependence [22]. Let us discuss this relevant issue in some more detail.

Schrödinger functional correlators
We now assume a finite volume with Dirichlet boundary conditions at x 0 = 0 and x 0 = T as explained in refs. [18,21,23]. These allow for the definition of gauge invariant boundary fields, for instance constructed from the boundary quark (and anti-quark) fields ζ (ζ) at x 0 = 0. In space we use periodic boundary conditions with a phase θ, Because of the zero momentum projection of the boundary fields, the correlator has a smooth behaviour at short distances (and small L) of the form in QCD with 3 colors and for massless quarks. There is no kinematical enhancement of a-effects. Together with the smooth background field introduced for the definition of the running coupling [23,24], this explains the very small lattice spacing effects in the running coupling and running operators mentioned before. We may also discuss the behaviour at large distances. There a saturation by few intermediate states in the spectral decomposition (see ref. [25] for a discussion of w, E) will give an accurate description of the correlation function. The energies E π n are the finite volume eigenvalues of the QCD Hamiltonian in the pion sector and E vac m are the energies of states with vacuum quantum numbers. Their dependence on L and the lattice spacing is suppressed. If the kinematics, given by θ, x 0 , L is such that effectively a few states with energies up to aE ≈ 1 contribute, the x 0 -dependence may again be strong and cutoff effects may be enhanced. However, the enhancement will not be as large as in Sect. 2.1 since there is no singular dependence at small x 0 . In Sect. 3 we will see quantitatively how the behaviour changes as L and x 0 are increased starting from L, x 0 ≪ 1/Λ.
We remark that the foregoing discussion is of course not in contradiction to the perturbative behaviour. In the perturbative region, x 0 is small and many states contribute significantly. But asymptotic freedom implies that their coefficients w nm are fine-tuned such as to produce the smooth behaviour of eq. (2.8). This is often called quark-hadron duality.

Energies and matrix elements
The most common application of Symanzik's effective theory is to energies and matrix elements, e.g. extracted from eq. (2.9) at large x 0 . The expansion of such observables, will be valid and accurate when the relevant energies, momenta and masses are small compared to a −1 . For the O(a)-improved theory, a first scaling test [15] showed that indeed the corrections to the continuum appear to be reasonably small when a ≤ 0.1 fm. We now proceed to discuss numerical results for f P in different regions of L (or g 2 (L)). Let us start from the short distance, weak coupling, region. In the nonperturbative computation of the scale-dependence of composite operators, such as P a , our collaboration chose Schrödinger functional boundary conditions with θ = 0.5, T = L and a vanishing background field [12]. From these simulations we have precise data for f P lat in a range of the Schrödinger functional coupling 0.98 ≤ḡ 2 (L) ≤ 5.5 [11,12]. In order to cancel the renormalization factor Z P , we normalize by f P lat at the midpoint, Atḡ 2 = 0, eq. (2.8) yields the slowly falling exponential plotted in Figure 1 as a continuous curve. The non-perturbative results atḡ 2 (L) ≈ 5.5 are next to this tree-level curve, even though in general perturbation theory is not very accurate at this coupling. For example the renormalization factor Z P is known to be a factor 2 away from one. We note in passing that the very close agreement with tree-level at x 0 /T ≥ 1/2 is somewhat accidental, as e.g. at a weak coupling ofḡ 2 (L/2) = 1.50 the curve is some 15 % away (dashed line). All in all, the smooth behaviour anticipated in our general discussion, is found to a remarkable degree in the whole range of L < ∼ 0.5 fm. Cutoff effects are very small. Even at the shortest distances, x 0 /a, the a 2 effects are only a few percent atḡ 2 (L) = 5.5. At weaker coupling,ḡ 2 (L/2) = 1.50, the a 2 effects are invisible within our precision of ≈ 2% (the data are not plotted in order not to clutter the graph). As explained in the previous section, this behaviour of the Schrödinger functional correlators goes hand in hand with small cutoff effects in the step scaling functions which describe for example the running of coupling and quark mass [12].
In physical units we have thus far investigated the region L = T < ∼ 0.5 fm. For L just somewhat larger, L = 0.8 fm, T = 9/4 L, and θ = 0, we have previously observed significant a 2 effects in our computation of Z A [10]. They were prominent in the statistically significant disconnected contributions to Z A , which can be shown to vanish in the continuum limit. As a consequence Z A determined with or without the disconnected diagrams differed by almost 15% at a ≈ 0.1 fm. Although the disconnected diagrams vanish quickly as a is reduced, an unpleasantly large ambiguity remained at typical values of a. For details we refer to the quoted reference.
We now turn to that same kinematics and subsequently to larger L. But first we point out that in this region the dependence of f P lat (x 0 )/f P lat (T /2) on θ is insignificant compared to the effects we will discuss. An explicit example is provided by comparing the dotted line of Figure 2 with the L/a = 8 data points. For our numerical demonstration we will thus freely use data at available values of θ.
We start with correlation functions f P lat at the same parameters as in [10], except for that we remain with θ = 0.5. The behaviour of f P lat , plotted in Figure 2, differs drastically from Figure 1 -see the tree-level curve as a reference. The non-perturbative correlators drop steeply (note that we use a logarithmic scale) and follow the characteristics of the described intermediate regime, where neither the smooth perturbative behaviour is realized, nor a single intermediate state is dominating (the latter would be seen as a straight line in the figure). For the coarsest lattice the logarithmic slope (effective mass) at x 0 = T /2 is as high as 0.8 in lattice units. In this situation, we find indeed very significant lattice spacing effects. It is then not surprising that there are also large a 2 effects in the form of the mentioned disconnected diagrams.
For an improved determination of Z A , we would like to choose a kinematical region, where (1) the a 2 effects are better suppressed and (2) the Schrödinger functional simulations can be done close to or in the chiral limit. The simulations can reach the chiral limit when the infrared cutoff on the spectrum of the Dirac operator is sufficiently large. Since in the Schrödinger functional the infrared cutoff is dominantly controlled by 1/T , we do not want to increase T significantly compared to the previous T ≈ 1.8 fm.
On the other hand, the energies E π n are expected to be decreasing with L; in fact in the small L limit they scale as L −1 . This leads us to consider a somewhat larger L, namely L ≈ 1.2 fm. Also in such a situation, namely with T = 3/2 L, we have simulation results [26], generated for the determination of improvement coefficients b m , b A −b P as well as the renormalization of the bare mass m q following ref. [27]. The circles in Figure 2 show f P lat in this kinematical situation 2 . The function has a much slower decay and indeed only moderate lattice artifacts. The energies of the states dominating the correlator appear to be significantly lower -at least the logarithmic slope never exceeds a value of 0.3 at the coarsest lattice spacing.
This kinematics is a good starting point for a redetermination of the renormalization factor Z A with reduced intrinsic a 2 ambiguities.

New determination of Z A
The renormalization condition of the axial current is obtained by considering an axial Ward identity in exactly the same way as it was done in ref. [10]. One starts from a Schrödinger functional correlation function of the axial current with two pseudo-scalar boundary operators, and performs an axial rotation of the variables in a region around the axial current [28]. Using isospin symmetry and PCAC one obtains for vanishing quark mass the Ward identity in terms of the correlation function The superscript I reminds us that the improved axial current is to be inserted. We refer to section 2 in ref. [10] for a derivation, the exact lattice implementation as well as a generalization to finite mass, which we use. Compared to that work we only change the kinematics. First of all the computation is now performed on a lattice of size L ≃ 1.2 fm instead of 0.8 fm, as motivated in the previous sections. Secondly, together with the sources O a , eq. (2.5), and O ′a , which correspond to ω 0 (x) = 1, we consider the following basis of wave-functions ω 1 (r) = e −r/a 0 , ω 2 (r) = r e −r/a 0 , ω 3 (r) = e −r/(2a 0 ) , in order to construct the external operators. Here, we keep the physical length scale a 0 fixed in units of L by choosing a 0 = L/6 and the (dimensionful) coefficients N i are set to normalize the wave function via a 3 x ω 2 i (x) = L 3 . The same set of interpolating fields, for the same L, has been used in ref. [29] to determine the improvement coefficient c A by requiring the quark mass derived from the PCAC Ward identity to stay the same as the external states are changed. We therefore choose the already determined [29] optimal wave-function i ω i , η (0) = (0.5172, 0.6023, 0.6081) . (4.8) It suppresses the contribution of the first excited state in the pseudoscalar channel to the correlation functions under consideration (see section 2 in ref. [29]). The final result will turn out to differ significantly from the determination in ref. [10] at the two largest couplings only. We therefore did not recompute Z A for small couplings and rather use the old estimates.

Results
We have two dynamical flavors of non-perturbatively improved Wilson quarks [30] and the plaquette gauge action. The improvement coefficients c sw , c A were set to their nonperturbative values [29,30]. We chose T = 3/2L with periodic boundary conditions (θ = 0) in space and vanishing background field. We used the HMC algorithm with two pseudo-fermion fields as proposed in refs. [1,31]. The particular implementation has been discussed and tested in refs. [32] and [33]. Following the last reference we chose a trajectory length of τ = 2 except for at β = 5.2 where we set τ = 1.
The normalization factor Z A is given by [10] The definition off I PA (2T /3, T /3, ω) and the PCAC quark mass m can be found in ref. [10]. After performing the Wick contractions one realizes that disconnected quark diagrams, where no propagator connects the two boundaries, contribute to f I XY (x 0 , y 0 , ω). These can be shown to vanish in the continuum massless limit as a consequence of the conservation of the axial current [10]. In an improved theory they therefore amount to O(a 2 ) effects on Z A . By dropping them one obtains an alternative definition of Z A , denoted Z con A . With the kinematics of ref. [10] the difference between Z A and Z con A for β = 6/g 2 0 < 5.5 was found to be rather large, though consistent with O(a 2 ) scaling. As it will become clear in the following this effect is very small in the computation presented here, which therefore significantly improves on the result in ref. [10].
In order to ensure a smooth dependence of Z A on the bare coupling g 2 0 and the correct scaling of discretization errors proportional to a 2 , we impose our normalization condition on a line of constant physics. This requires keeping all length scales fixed as g 2 0 is varied. The lattice size L is set to approximately 1.8L * with L * given by the conditionḡ 2 (L * ) = 5.5, whereḡ 2 (L) is the Schrödinger functional coupling defined in refs. [11,24]. The relation between L * a and g 0 could be taken from ref. [34]. The bare parameters of our simulations and the results for Z A are collected in Table 1.
Due to algorithmic instabilities caused by the appearance of very small,   unphysical, eigenvalues in the spectrum of the Wilson-Dirac (SF) operator [3,35], at the coarsest lattice spacing we could simulate down to bare quark masses of about a m ≈ 0.01 only. Our estimates for the massless limits at β = 5.2 are just weighted averages of the numbers at the two lightest quark masses. The errors in Table 1 include a systematic uncertainty for possible O(am) contaminations given by the difference between the determination at the heaviest mass and the described estimate. This uncertainty is added linearly to the statistical error. A similar O(am) uncertainty would be obtained by comparing to a linear fit in am with all three points. Figure 3 illustrates the procedure for Z con A and ω = ω π (0) . It is clear from the plot that the dependence of Z A on the quark mass is rather mild as expected for the "massive"definition of the renormalization constant [10] used here. Similar remarks obviously apply to the case ω = ω 0 , as Table 1 shows that the two wave-functions give consistent results for all values of κ. At the other values of the bare coupling we simulated at very small quark masses. Given also the flat dependence observed at β = 5.2 we did not need to estimate an effect of order am. As anticipated, already at β = 5.7 the present result nicely agrees with the numbers in ref. [10]. Conversely, at β = 5.2 the new determination of Z A is about 8% larger than the old one. These comparisons are summarized in Figure 4.
In view of the above discussion we smoothly interpolate the data for Z con A and ω = ω π (0) in the region 5.2 ≤ β ≤ 5.7 and those in ref. [10] for 7.2 ≤ β ≤ 9.6 (filled circles in Figure 4) by the formula where the coefficient of the term linear in g 2 0 is fixed by 1-loop perturbation theory [36] and the last two coefficients are the result of a fit. We ascribe an absolute error to Z A , with the old ones for ZA and Z con A in ref. [10]. The solid curve represents the interpolation formula in eq. (4.10).

Other renormalization conditions of the axial current
In order to get a further impression about residual cutoff effects, we also studied an alternative definition of Z A . It is obtained in the same framework by replacing a light quark with a static one in the external operators O and O ′ . By remaining with the flat wave-function and by denoting the static quark field on the boundary x 0 = 0 by ζ h we write where i = 1, 2 is a flavour index. The flavour contractions in f I XY are changed correspondingly. The correlator f hl,I AA for example is written and the massless normalization condition becomes with (4.14) As the fields X and Y in the correlator do not contain static fields, it is clear that disconnected diagrams cannot appear and the static quark propagates from one boundary to the other. Static quarks are discretized through the HYP1 and HYP2 static-quark actions [37], which have a relatively good signal to noise ratio in static-light correlation functions at large time separations.  The check was performed at the largest lattice spacing where ambiguities are expected to be most pronounced. Simulation parameters and results are collected in Table 2. The results at L ≃ 1.2 fm agree with those in the previous section, indicating again small overall cutoff effects. The values at L ≃ 0.8 fm instead again suggest that the determination in ref. [10] suffers from large a 2 effects, as both Z A and Z con A there differ significantly (and especially for small quark masses) from the numbers in Table 2.

4.3
Renormalization of the vector current We also recomputed the renormalization constant Z V of the vector current in the new kinematics. Since it changes by less than 2% compared to ref. [10] at the largest lattice spacing, there is no reason to publish a new determination.

Conclusions
We have discussed possible kinematical enhancements of cutoff effects in lattice gauge theory determinations of QCD correlation functions. In this respect, the typical SF correlation functions improve significantly over those of composite local fields with periodic boundary conditions (or large volume). The origin of this difference is simply the difference in mass dimensions of the Schrödinger functional boundary fields compared to the usual local composite fields. Since the latter is at least three, time-slice correlators diverge at least as x −3 0 at short distances. Such a steep behavior is difficult to approximate in a discretized theory.
In the Schrödinger functional we identified a parametrically much weaker, but still relevant, kinematical enhancement of cutoff effects. It appears in the transition region between approximately perturbative behaviour and dominantly non-perturbative one. Numerically we find that it appears e.g. for a T × L 3 geometry with L ≈ 0.8 fm , T ≈ 1.8 fm. In this kinematical situation a few hadronic intermediate states are relevant for the correlators and can produce a relatively steep decay of correlation functions even at a small quark mass.
This means that the transition region from approximately perturbative to strongly non-perturbative is a relatively difficult one for numerical simulations. Discretization errors have to be investigated carefully. Fortunately we also saw that this region is rather narrow. With the step scaling method [19,20,22] it is typically bridged by one step.
Avoiding the difficult region in the renormalization condition for the relativistic axial current, we have finally presented a significant improvement of the previous determination of Z A [10]. The difference to the old one is of order a 2 , but it is up to 10 % at the largest lattice spacing considered. In our new determination we find that disconnected contributions, which can be shown to vanish in the continuum limit, are very small -in contrast to the previous kinematical setup [10]. We also see a nice agreement with a renormalization condition where a static quark is present as a spectator in the Ward identity. The previous determination of Z V is confirmed.