High precision renormalization of the flavour non-singlet Noether currents in lattice QCD with Wilson quarks

We determine the non-perturbatively renormalized axial current for O($a$) improved lattice QCD with Wilson quarks. Our strategy is based on the chirally rotated Schr\"odinger functional and can be generalized to other finite (ratios of) renormalization constants which are traditionally obtained by imposing continuum chiral Ward identities as normalization conditions. Compared to the latter we achieve an error reduction up to one order of magnitude. Our results have already enabled the setting of the scale for the $N_{\rm f}=2+1$ CLS ensembles [1] and are thus an essential ingredient for the recent $\alpha_s$ determination by the ALPHA collaboration [2]. In this paper we shortly review the strategy and present our results for both $N_{\rm f}=2$ and $N_{\rm f}=3$ lattice QCD, where we match the $\beta$-values of the CLS gauge configurations. In addition to the axial current renormalization, we also present precise results for the renormalized local vector current.


Introduction
Lattice regularizations with Wilson type fermions [3] are widely used in current lattice QCD simulations [4][5][6][7][8][9][10]. The ultra-locality of the action enables numerical efficiency and thus access to a wide range of lattice spacings and spatial volumes. Furthermore, Wilson fermions maintain the full flavour symmetry of the continuum action, as well as the discrete symmetries such as parity, charge conjugation and time reversal. Unitarity is either realized exactly, or, in the case of Symanzik-improved actions, approximately up to cutoff effects which vanish in the continuum limit.
The price to pay for these advantages consists in the explicit breaking of all chiral symmetries by the Wilson term in the action. Well-known consequences include the additive renormalization of quark masses, the mixing under renormalization of composite operators in different chiral multiplets and discretization effects linear in a, the lattice spacing. Furthermore, the Noether currents of chiral symmetry are no longer protected against renormalization.
The matrix elements of the axial Noether currents between pion or kaon states and the vacuum, parametrized by the decay constants f π,K , e.g.
can be related to the measured life times of pions and kaons. The decay constants are finite in the chiral limit, can be precisely measured in numerical simulations and are ideally suited to set the scale in physical units. In order to achieve this with Wilson quarks one needs to determine the correctly renormalized axial currents, (with flavour indices f 1,2 = u, d, s), which are to be inserted into the matrix elements. Of course it is desirable that the error of the matrix elements is not dominated by the uncertainty of the current normalization constant. Over the last 30 years many efforts have been made to control the consequences of explicit chiral symmetry breaking with Wilson quarks. The main strategy consists in imposing continuum chiral symmetry relations as normalization conditions at finite lattice spacing [11,12]. This is usually done using chiral Ward identities, which follow from an infinitesimal chiral change of variables in the QCD path integral. An example is the PCAC relation which determines the additive quark mass renormalization constant, as the "critical value" of the bare mass parameter, where the axial current is conserved. The fact that chiral symmetry is fully recovered only in the continuum limit implies that the choice of normalization condition matters at the cutoff level; at a fixed value of the lattice spacing the numerical results may occasionally differ substantially between any two such choices. Rather than interpreting this scatter as a systematic error, the modern approach consists in choosing a particular normalization condition and in fixing all dimensionful parameters (such as momenta or distances or background fields) in terms of a physical scale. This defines a so-called "line of constant physics" (LCP), along which the continuum limit is taken. As the lattice spacing a (or, equivalently, the bare coupling, g 2 0 = 6/β), is varied, this defines a well-defined function Z A = Z A (β). Obviously, another choice for the LCP will result in a different function Z A (β). However, their difference will be, within errors, a smooth function of β which vanishes asymptotically ∝ a or ∝ a 2 if O(a) improvement is implemented. Hence, following a LCP ensures that cutoff effects are smooth functions of β and the choice of LCP becomes irrelevant in the continuum limit. Adopting this viewpoint, the relevant systematic error is therefore determined by the precision to which a chosen LCP can be followed.
In this paper we apply a recently developed method to lattice QCD with N f = 2 and N f = 3 flavours, matching the lattice actions chosen by the CLS initiative [8,10]. Our method is based on the chirally rotated Schödinger functional (χSF) [13,14]. The theoretical foundation of this framework has been explained in [14] and it has passed a number of perturbative and non-perturbative tests [15][16][17][18][19][20]. In contrast to the Ward identity method the axial current renormalization conditions follow from a finite chiral rotation in the massless QCD path integral with Schrödinger functional (SF) boundary conditions. The renormalization constants are then obtained from ratios of simple 2-point functions. For the axial current, this represents a significant advantage over the Ward identity method [12,21,22] which involves 3-and 4-point functions. Hence, we observe a dramatic improvement in the attainable statistical precision for Z A and some care is required to ensure that systematic errors are under control at a similar level of precision. We also discuss the normalization procedure for the local vector current. While flavour symmetry remains unbroken on the lattice with (mass-degenerate) Wilson quarks, the corresponding Noether current lives on neighbouring lattice points connected by a gauge link, so that the use of the local vector current is often more practical.
This paper is organized as follows: after a short reminder of the χSF correlation functions in the continuum and the normalization conditions derived from them in sect. 2, we define in sect. 3 a couple of different LCPs which we have followed. We then present the Z A and Z V determinations for lattice QCD with N f = 2 and N f = 3 quark flavours in sects. 4 and 5, respectively, together with various tests we have carried out. Sect. 6 contains a summary of the main results of this work and some concluding remarks. Finally, the paper ends with two technical appendices: appendix A collects the parameters and results of the simulations, and appendix B provides a detailed discussion on the systematic error estimates for our determinations.
The main results for N f = 2 are collected in table 4, while those for N f = 3 are given in tables 6 and 7. These results can be directly applied to data obtained from the CLS 2-and 3-flavour configurations, respectively [8,10]. The N f = 3 results have, in fact, already been used, and enabled the precision CLS scale setting in ref. [1] and the accurate quark-mass renormalization of ref. [23].

The Schrödinger functional and chiral field rotations
We start by considering massless two-flavour continuum QCD. The Euclidean space-time is taken to be a hyper-cylinder of volume L 4 with Schrödinger functional boundary conditions [24,25]. In particular, in the Euclidean time direction, the quark and anti-quark fields satisfy, and similarly at time x 0 = L with the change P ± → P ∓ . The SU(2)×SU(2) chiral and flavour symmetry leads to conserved iso-vector Noether currents, given by A a µ (x) = ψ(x)γ µ γ 5 τ a 2 ψ(x), V a µ (x) = ψ(x)γ µ τ a 2 ψ(x), (2.2) with Pauli matrices τ a and isospin index a = 1, 2, 3. SF correlation functions of these currents with iso-vector boundary sources O a 5 and O a k have been defined in [26,27] and are given by Passing from the isospin notation to fields with definite flavour assignments, and similarly for the boundary sources, the correlation functions for isospin indices a = 1, 2, can be written in terms of the flavour off-diagonal fields, For the flavour diagonal fields in the isospin a = 3 components, e.g.
one may use flavour symmetry to write and analogously for k V . Note that the additional up-type flavour u is merely a notational device to indicate the fermionic contractions taken into account when applying Wick's theorem. Indeed, the sum of all the disconnected contributions for the flavour diagonal a = 3 components of SF correlation functions cancels exactly due to flavour symmetry. We now apply a flavour diagonal chiral rotation to the fields, Choosing the rotation angle α = π/2 then leads to the chirally rotated SF boundary conditions,Q with projectorsQ ± = 1 2 (1±iγ 0 γ 5 τ 3 ). Analogous boundary conditions with reverted projectors are obtained at x 0 = L. Applying the same chiral field rotation to the axial currents, one obtains either a vector current or remains with an axial current, depending on the flavour assignments. If the chiral rotation of the field variables is performed as a change of variables in the functional integral, one arrives at the formal continuum identities where the gand l-functions are defined with χSF boundary conditions, eqs. (2.9), for instance Here, the boundary operators Q f 1 f 2 5 denote the chirally rotated versions of their SF counterparts, O f 1 f 2

5
. For the complete expressions and further details we refer to ref. [20]. Regarding the case of QCD with N f = 3 quark flavours we note that the very same steps can be taken provided the massless third quark does not take part in the chiral rotation and thus remains with standard SF boundary conditions [14]. Correlation functions are then considered for the doublet fields only, i.e. the third quark never appears as a valence quark.

Renormalization conditions
In the lattice regularized theory with Wilson type quarks, relations such as (2.11) can only be expected to hold after renormalization and up to cutoff effects. One first has to ensure that massless QCD with χSF boundary conditions has been correctly regularized. This is achieved by tuning the bare mass parameter m 0 to its critical value, m cr , where the axial current is conserved, and by tuning a boundary counterterm coefficient z f such that physical parity is restored (cf. [20] for more details). In terms of the bare χSF correlation functions one may choose the two conditions, (with∂ 0 the symmetric lattice derivative). The division by the pseudo-scalar correlation function g ud P is not really necessary, however it is done for convenience, as it gives rise to the definition of a (bare) PCAC quark mass m. Solutions to these equations define m cr and z * f as functions of the bare coupling g 0 , and the lattice size, L/a. Once the lattice regularization is correctly implemented, we expect e.g.
where Z ζ renormalizes a boundary quark or anti-quark field [14,26,28] and Z A,V are the current normalization constants of interest. Requiring such identities to hold exactly at finite lattice spacing thus fixes the relative normalization of axial and vector current. Replacing the latter by the exactly conserved lattice vector current V µ (x) (cf. ref. [20]), for which Z V = 1, one may obtain Z A from either one of the ratios Assuming that the parameters x 0 (here set to L/2), the boundary angle θ [29], the background gauge field [24], and the precise definition for the zero mass point (2.13) are fixed, we define, on an (L/a) 4 lattice and for a given bare coupling g 2 0 = 6/β, Finally the choice of a line of constant physics (cf. section 3) defines a smooth function (L/a)(β) such that the normalization constants become functions of β alone, with the difference between any two definitions vanishing smoothly with a rate ∝ a 2 . We also comment on the appearance of a second up-type flavour u in (2.15). When applying the chiral rotation (2.8) to the diagonal components of f A , the disconnected diagrams are mapped to disconnected diagrams on the χSF side which can be shown to add up to a pure cutoff effect. Their omission is thus perfectly legitimate, even if the formulation of the renormalization conditions then has an element of partial quenching to it. The situation is comparable with the Ward identity method in two-flavour QCD [12,21], where a fictitious s-quark can be introduced to eliminate the disconnected diagrams.
Even though there exists a conserved vector current, in practice the local current is often used and then requires renormalization, too. Its renormalization constant can be obtained from, The same remarks as for the axial current normalization apply here, and with definite choices for all parameters we set, As in the case of the axial current normalization conditions, only 2-point functions are required, which connect the boundary quark bilinear sources with the currents in the bulk. This is a major advantage over the Ward identity method [12,21] where 3-and 4-point functions are required. Hence, one expects better statistical precision from the simpler 2-point functions, and this will be confirmed below. Furthermore, as discussed in [20], the cutoff effects in the ratios are O(a 2 ), due to the mechanism of automatic O(a) improvement, even if the PCAC mass and the axial current are not O(a) improved by the counterterm ∝ c A [26], or if the vector currents are not improved by the corresponding counterterms ∝ c V , c V [27,30]. Finally, we emphasize that similar renormalization conditions can be devised for other finite renormalization constants. An interesting example is the ratio Z P /Z S , where Z P and Z S are the pseudo-scalar and scalar renormalization constant, respectively. We refer the reader to ref. [20] for more details.

General considerations
A line of constant physics requires to specify a physical (length) scale r which is kept fixed as the continuum limit is taken. A typical choice would be the pion decay constant, r = 1/f π , either at the physical quark masses or in the chiral limit. Once calculated for a range of lattice spacings, this scale defines a function (r/a)(β) of the bare coupling β = 6/g 2 0 which fixes the lattice spacing a in units of the chosen physical scale. Choosing the spatial lattice extent L/a, at a given beta, such that (L/a)(β) (r/a)(β) = L/r = C r (3.1) (with a numerical constant C r ) then fixes the spatial size of the finite volume system in units of r. In practice we will choose C r such that the physical size of L will be somewhat larger than half a femto metre. Note that this equation can be read in two ways: first, if one fixes C r and then chooses a set of β-values for which r/a is known, one obtains a corresponding set of values (L/a)(β), which will not necessarily be integers. To evaluate the normalization constants at these non-integer lattice sizes then requires some interpolation of results from neighbouring integer L/a-values at the same β. Alternatively, one could choose a set of integer L/a-values such that a choice for C r will imply a set of β-values.
In general this means that the data for r/a may have to be interpolated in β. We will here choose the first option, with the set of β-values taken over from the large volume simulations by the CLS project [8,10].
Having set the scale one needs to ensure the correlation functions are calculated in the desired situation of massless QCD and for the chosen chirally rotated boundary conditions at α = π/2. This means one needs to tune the bare quark mass am 0 and z f as functions of β. We will discuss this in more detail below. Finally, the correlation functions depend on kinematical parameters, such as x 0 or background field parameters such as θ. We have already set x 0 = L/2 in eqs. (2.15,2.17) and we choose θ = 0 and work with vanishing SU(3) background field.
With these parameter choices we will have, for a given r and C r in eq. (3.1), two definitions each for Z A and Z V , namely either based on the gor the l-ratios. We then expect e.g. that where the a 2 -effects are now expected to be smooth functions of the bare coupling.

Perturbative subtraction of cutoff effects
A possible refinement consists in using perturbation theory to reduce the cutoff effects perturbatively. This requires to compute the R-ratios (2.15,2.17) perturbatively, with the exact same parameter choices as in the numerical simulations. We have performed this calculation to 1-loop order, and for the chosen parameters we always find R g,l(0) A,V (a/L) = 1 exactly. We may then define a 1-loop correction factor, and results for the coefficients R g,l(1) A,V are collected in table 1, for the relevant lattice resolutions L/a and the two lattice gauge actions used by CLS. Note that the 1-loop results are N f -independent and are thus obtained along the lines of ref. [20], the only difference being the form of the free gluon propagator in the case of the Lüscher-Weisz gauge action [44]. As an aside we note that our results converge to the known 1-loop results Z (1) A,V for an infinitely extended lattice [31][32][33], i.e. for a/L = 0. We also observe that the 1-loop cutoff effects for the l-definitions are generally much smaller than for the g-definitions.
Wilson gauge action  For given L/a and β, the perturbatively improved current normalization constants are now defined by Z g,l A,V, sub (β, L/a) = r g,l A,V (β, L/a) × Z g,l A,V (β, L/a), (3.6) and, by construction, the O(a 2 ) cutoff effects are subtracted to O(g 2 0 ), reducing them to O(a 2 g 4 0 ). The subtracted data for the Z-factors are then treated as before: a choice of a line of constant physics implies a set of βand corresponding L/a-values to which the data must be interpolated. We will see evidence for the effectiveness of this perturbative subtraction of cutoff effects in Sect. 4 and 5.
3.3 Choices of LCP for N f = 2 and N f = 3 In order to fix the physical scale r, we choose either the kaon decay constant r = 1/f K (N f = 2), or the gradient flow scale r = √ 8t 0 (N f = 3) [34]. 1 In order to fix the respective constants C r we proceed as follows. Given the set of values β i for i = 1, 2, . . . (taken from CLS), we choose as a reference value β ref either the largest or the smallest of the set. Choosing an integer lattice size L/a at the reference point β ref now fixes C r through Having set the scale in this way, the L/a-values at the remaining β i follow from eq. (3.1). For all our choices the physical size of our space-time extent will be L ≈ 0.6 − 0.7 fm. As mentioned before, except at the chosen reference value for β this requires interpolations of simulation results at integer L/a and our current simulation code, which is based on the openQCD package [35,36], requires that L/a is also even.

Topology freezing
Numerical simulations of the SF by means of standard Monte Carlo algorithms are known to suffer from the topology freezing problem (see e.g. ref. [37] for a discussion). A possible solution is to follow the proposal of ref. [37] and simulate the theory with open-SF boundary conditions. However, if for the given choice of parameters the problem is "mild", one can circumvent the issue in a straightforward manner by simply imposing the renormalization conditions (2.18) and (2.13) within the trivial topological sector [38,39]. In a continuum notation, the correlation functions entering these definitions are modified as follows, , (3.8) and analogously in all other cases. 2 Here, the Kronecker δ in the functional integral selects the gauge field configurations with topological charge Q = 0. Since relations based on chiral flavour symmetries should hold separately in each topological charge sector, this restriction to the trivial sector is a legitimate modification of the current renormalization conditions. It provides a viable solution to the algorithmic problem of topology freezing in cases where this problem becomes marginally relevant; this means when the fraction of topologically non-trivial gauge field configurations in the relevant ensembles is not too large. For our choices of parameters, the percentage of gauge field configurations with Q = 0 is generally below 10%, and reaches approximately 30% only in a couple of cases (cf. tables 8 and 9). On the lattice the topological charge is not unambiguously defined. We follow refs. [38,39] and define the trivial topological sector as the set of gauge field configurations for which |Q| < 0.5, where Q is discretized in terms of the Wilson flow and the clover definition of the field strength tensor [34]. The flow time t is then kept fixed in physical units by requiring √ 8t = 0.6 × L.

On the tuning of am 0 and z f
The current normalization conditions require the χSF correlation functions at zero quark mass and with a chiral twist angle of π/2. In practice this is achieved by the simultaneous tuning of m 0 and z f such that eqs. (2.13) are satisfied. In general a 2-parameter tuning can be quite involved. However, here the non-perturbative O(a) improvement of the action implies that the O(a) uncertainty of the zero mass point is very much reduced. Since a change in z f merely re-defines the matrix element used to define the PCAC mass, a variation of z f is expected to induce a small variation of m within this O(a) uncertainty. The latter could in principle be reduced to O(a 2 ) by including the c A -counterterm to the axial current, but this will not be pursued here. Another important observation is that, once m 0 and z f are within O(a) of their target values, the sensitivity of the PCAC mass to a variation of z f is reduced to order a 2 (cf. appendix B, discussion after eq. (B.8)). One is therefore led to conclude that the PCAC mass m is to a good approximation independent of z f , and the tuning of m 0 and z f thus becomes straightforward; given a reasonable guess for z f , one can first tune m 0 , and then turn to z f . As an illustration of this situation we discuss the N f = 2 case for L/a = 8, β = 5.3. For the tuning we considered 3 values of κ = 1/(2am 0 + 8) and 4 values of z f . We then generated around 2000 gauge field configurations separated by 10 MDUs for each of the 12 ensembles, and measured the relevant correlation functions. Figure 1 collects the results for the PCAC mass as a function of the bare quark mass, for the 4 different values of z f . Within statistical errors, the PCAC mass depends linearly on m 0 and is essentially independent of z f . A linear fit of m vs. m 0 yields an estimate of m 0 = m cr (g 0 , L/a) for which m vanishes: these are collected in table 2. The results are perfectly compatible with each other, and we take as our estimate for m cr the result of a weighted average of these four.
Once the critical bare mass is fixed, a smooth interpolation of g ud A (L/2) in m 0 gives the results shown in figure 2. Over the chosen range, g ud A (L/2) so interpolated is perfectly linear in z f , and it is thus straightforward to determine the point z * f where g ud A (L/2) vanishes i.e. z * f = 1.2877 (5) in this example.   The estimated values of am cr and z * f determined in this way turn out to be quite accurate in practice, cf. table 8. 3 We remark that results for m cr could also be taken from a different source, for instance from standard SF simulations. In this case only z f needs to be tuned. The differences to the above procedure would be O(a) both in m cr and in z * f which, by the mechanism of automatic O(a) improvement, induce O(a 2 ) differences in observables such as the current normalization constants [14,20]. One also expects that a precise tuning of m 0 is less crucial in the χSF than in the SF; the quark mass dependence of physical observables around the chiral limit is quadratic rather than linear [40].

Sources of uncertainties
Besides statistical errors directly affecting the estimators for the current normalization constants, the other source of uncertainty originates from the precision to which a line of constant physics can be followed. In principle also this latter effect is of a statistical nature, however, some elements of modelling or estimates may be involved when propagating these errors to the normalization constants, so that it is partly justified labelling these effects as systematic.
Our procedure consists of the following steps: 1. The LCP together with the set of values β i translates to target values (L/a)(β i ).
At each β i we choose lattices with even L/a straddling the target values. We here anticipate that with our choices of LCPs the required lattice sizes are in the range L/a = 8 to L/a = 16. Note that all target values (L/a)(β i ) come with statistical errors except for β = β ref , where, by definition, L/a is given as an (even) integer.
2. For given β and L/a we determine the solutions am 0 = am cr and z f = z * f of eqs. (2.13). In order to find their statistical errors which follow from the statistical uncertainties on m and g ud A (L/2), we use estimates for the relevant derivatives, 3. We then determine the induced error on the Z-factors by estimating their derivatives with respect to the bare parameters, It turns out that the derivatives (3.9) and (3.10) scale quite well with lattice size and lattice spacing, so that it is unnecessary to evaluate them for all parameter choices. Some cross checks are sufficient. The errors coming from the uncertainties in m and z f are then combined in quadrature and added, again in quadrature, to the statistical error.
estimate for the derivative is used to assign a systematic error due to the difference ∆(L/a) ≡ L/a − (L/a)(β i ), also taking into account the statistical uncertainty on (L/a)(β i ). The resulting systematic error is again added in quadrature.
We emphasize that all systematic effects become essentially statistical errors provided enough data is produced to estimate the derivatives required to propagate the errors to the normalization constants. In the following two sections we will present the lattice set-up and results for N f = 2 and N f = 3 lattice QCD. We will also come back to some of the above points. 4 Numerical results for N f = 2 flavours

Lattice set-up and parameter choices
The CLS large volume simulations of 2-flavour QCD [8] were performed using non-perturbatively O(a) improved Wilson quarks and the Wilson gauge action. The matching to CLS data via the bare coupling requires that we use the same action in the χSF. As for the details of the action near the time boundaries we refer to ref. [20]. In particular the counterterm coefficients c t (g 0 ) and d s (g 0 ) were set to their perturbative one-loop values using the results of that reference. In general, the incomplete cancellation of boundary O(a) artifacts implies some remnant O(a) effects in observables. However, for the estimators of the current normalization constants, eqs. (2.15,2.17) it can be shown that such O(a) effects only cause O(a 2 ) differences [20]. The CLS simulations were carried out for 3 values of the lattice spacing [8], corresponding to the β-values 5.2, 5.3 and 5.5. For future applications we have added a finer lattice spacing corresponding to β = 5.7. We choose the smallest CLS-value β = 5.2 as reference value and set L/a = 8 at β = 5.2, (4.1) to define the starting point for the line of constant physics. We then fix the space-time volume of the χSF simulations in terms of the kaon decay constant, f K , evaluated at physical quark masses. Taking af K from table 3 at β = 5.2 yields and corresponds to L ≈ 0.6 fm. Imposing this condition at the other β-values then leads to the (non-integer) L/a-values given in table 3. The quoted errors are a combination of statistical uncertainties, propagated from eq. (4.2) and the error on af K at the given β's. While the first 3 results for af K in table 3 have been directly measured [8] we have estimated af K at the fourth value, β = 5.7, as follows: with af K at β = 5.5 taken as starting point we used the three-loop β-function for the bare coupling [41], in order to determine the ratio of lattice spacings. The error is obtained by summing (in quadrature) the statistical error propagated from the result at β = 5.5, and a systematic error due to the use of perturbation theory. The latter is estimated as the difference between the non-perturbative result for af K at β = 5.5, and the same perturbative procedure, applied

Results and error budget
In table 4 we collect the results for Z A,V , both g and l definitions, at the four values of the lattice spacing. The statistics range from 1, 800 to 12, 000 measurements depending on the ensemble, cf. table 8. The quoted uncertainties combine the statistical and systematic errors. The statistical errors are at the level of 0.1 − 0.4 , depending on the Z-factor and ensemble considered. Hence a significant contribution to the error comes from systematic uncertainties. As discussed in section 3.6, systematic errors result from uncertainties or deviations in following a chosen LCP, which correspond with statistical errors and deviations from zero in m and g ud A (L/a), as well as uncertainties in the target lattice extent L/a and systematic errors arising from inter-or extrapolations from the simulated lattices sizes, if applicable. Tables 3 and 8 contain the relevant information for the case N f = 2. The propagation of these uncertainties to the Z-factors is then performed following the steps outlined in Sect. 3.6. We have carried out some additional simulations to estimate the derivatives in eqs. (3.9,3.10), and some perturbative calculation to check the expected scaling of the derivatives with the lattice size. We delegate a detailed discussion to appendix B. Here we just note that with our statistics and our rather conservative approach, the propagated uncertainties are typically larger than the statistical errors for the R-estimators eqs. (2.15,2.17) (cf. table 10 and 11).

Effect of perturbative one-loop improvement
As discussed in section 3.2, we have also computed the relevant χSF correlation functions in perturbation theory to order g 2 0 = 6/β. Besides consistency checks and qualitative insight the main application consists in the perturbative subtraction of cutoff effects from the data. Note that this requires to emulate the non-perturbative procedure in all details, in particular the determination of am cr and z * f according to eqs. (2.13). The lower part of table 4 contains the results for Z A,V after perturbative improvement. Comparing with the unimproved results in the upper part of table 4, one can see that the g-definitions are more affected, and are brought closer to the corresponding l-definitions by the perturbative improvement (cf. also figure 5). In any case, the perturbative corrections are at the level of 2 per cent at most. In conclusion, our final results for Z A,V , either with or without perturbative improvement, turn out to be very precise and improve significantly on the standard SF determination based on chiral Ward identities [8,21,42]. This is particularly true for the case of Z A , which can be appreciated in figure 3 where the determinations of table 4 are compared with those of refs. [8,42]. In figure 4 we show instead a comparison for the case of Z V , as obtained from the χSF, cf. table 4, and from the standard SF (cf. ref. [21]). We note that a relevant contribution to the error of our results comes from propagating the uncertainties associated with maintaining the condition (4.2) i.e. keeping L constant (cf. table 10 and 11). We anticipate that due to the much more accurate knowledge of the LCP in terms of t 0 (cf. table 5), and by using interpolations in L/a at all relevant β values, this source of error will be essentially eliminated in the case of N f = 3 (cf. section 5).

Universality and automatic O(a) improvement
The χSF determinations (2.15) and (2.17) are expected to be automatically O(a) improved once the bare parameters m 0 and z f are properly tuned (cf. section 2.2). This means that neither bulk nor boundary O(a) counterterms are necessary to cancel O(a) discretization errors in these quantities. This was confirmed to one-loop order in perturbation theory [20]   The individual SF points are taken from refs. [8,42], and are slightly displaced on the x-axis for better clarity. The solid black line corresponds to the SF results from the fit formula of ref. [8], and the dashed lines delimit the 1σ region of the fit. Note that the SF fit formula is obtained by considering additional points with g 2 0 < 1, here not shown, and by enforcing the perturbative 1-loop behaviour for g 2 0 → 0 (see ref. [8] for the details). The individual SF points are taken from refs. [21], and are slightly displaced on the x-axis for better clarity. The solid black line corresponds to the SF results from the fit formula of ref. [21], and the dashed lines delimit the 1σ region of the fit. Note that the SF fit formula is obtained by considering additional points with g 2 0 < 1, here not shown, and by enforcing the perturbative 1-loop behaviour for g 2 0 → 0 (see ref. [21] for the details).  and should hold generally. To this end we now look at the ratios between Z-factors coming from the gand l-definitions. The expectation that these ratios converge to 1 with O(a 2 ) corrections is indeed very well borne out by the data, cf. figure 5, where we also include fits to this expected behaviour. We emphasize that this is a non-trivial result: even though the bulk action is improved to match the CLS set-up, we did not O(a) improve the currents entering the definitions (2.15,2.17) and (2.13). This result thus confirms automatic O(a) improvement at the non-perturbative level, and, indirectly, the universality relations between the χSF and SF formulations. A direct way to test universality between the χSF and SF formulations would be simply to study the continuum scaling of ratios of Zfactors as obtained from one and the other formulation. Provided the SF determinations are properly improved, these should also approach 1 in the continuum limit with O(a 2 ) corrections. The large errors on the SF determinations do not allow us for a precise test of this expectation. However, the results in figure 3 and 4 clearly show that our determinations are in fact compatible with the SF ones within errors.  [1,2,10]. For completeness we note that CLS has also tried to simulate at a coarser lattice spacing corresponding to β = 3.3. However, these ensembles have been discarded for the scale determination in [1] due to very large cutoff effects observed e.g. in t 0 [10]. For this reason we will not consider this β-value in our study, however, we mention that it was adopted as starting point for the Ward identity determination of Z A in ref. [22]. Given the relatively large set of lattice spacings we here consider two different LCPs, with slightly different physical extent, L 1 and L 2 , which we define through the gradient flow time t 0 [34]. The associated length scale r = √ 8t 0 can be interpreted as a smoothing radius, and has been very precisely determined for the CLS β-values ≥ 3.4 in [1]. Using this scale we impose the conditions (16)  respectively. Using the result for t 0 in physical units [1], eqs. (5.1) translate to L 1 ≈ 0.7 fm and L 2 ≈ 0.6 fm.   In table 5 we collect the relevant β values of the CLS simulations and the corresponding results for t 0 /a 2 [2]. The latter are evaluated for equal up-, down-, and strange-quark masses, which are close to the physical average quark mass (see refs. [1,2]). Table 5 also gives the lattice sizes (L 1,2 /a)(β) which satisfy the conditions (5.1). Compared to the N f = 2 case (cf. table 3), it is obvious that these N f = 3 LCPs are much more accurately determined. In order to exploit this higher precision, we performed simulations for several L/a-values at each β (cf. table 5). This allowed us to accurately interpolate the Z-factors to the target values (see appendix B.5 for more details). Table 9 contains a summary of all simulations performed with the corresponding parameters. Due to both technical and historical reasons, we do not use the finest lattice spacing for the LCP defined in terms of L 1 . Following this LCP up to β = 3.85 would have required simulating lattices with L/a = 18, 20, which are particularly inconvenient to parellelize with our current simulation program. Note also that CLS simulations at β = 3.85 are ongoing and currently limited to a single ensemble, so that the LCP with L 1 may remain useful for a while. More importantly, however, the comparison between both LCPs allows us to perform additional tests on our results (cf. section 5.3).
The lattice action we employ for the finite volume simulations matches the CLS action in the bulk, i.e. the Lüscher-Weisz tree-level improved gauge action and 3 flavours of nonperturbatively improved Wilson quarks [43]. Close to the time boundaries of the lattice there is some freedom regarding the implementation of Schrödinger functional boundary conditions. For the gauge fields we choose option B of ref. [44]; we refer the reader to this reference for the details. Regarding the fermions, two quark flavours satisfy χSF boundary conditions (option τ = 1 of [14]), while the third one obeys the standard SF boundary conditions [25]. In general, such a mixed set-up increases the number of O(a) improvement coefficients which need to be tuned in order to eliminate O(a) discretization errors from the time boundaries. As in the N f = 2 case, however, one can show that the corresponding counterterms affect the renormalization constants Z A,V only at O(a 2 ). For definiteness we have used the one-loop estimate c t = 1 + g 2 0 c (1) t , where the one-loop coefficient decomposes as follows, c The pure gauge contribution taken from ref. [45], the fermionic χSF contribution from ref. [20] and the SF contribution from ref. [29]. 4 Furthermore, we use the tree-level values d s = 1/2 [20] andc t = 1 [26].

Results and error budget
In table 6 and 7 we collect the results for Z A,V , corresponding to the L 1 -and L 2 -LCP, respectively. The statistics we accumulated for the different ensembles ranges between 3,200 and 31,000 measurements, with exact numbers given in table 9. The corresponding statistical precision on the Z-factors is between 0.1 − 0.55 , depending on the exact quantity and ensemble. The errors quoted in the tables then combine the statistical errors with the systematic errors originating from the uncertainties on the LPCs. Like in the  N f = 2 case, the high statistical precision requires a careful assessment of the systematic errors in order to arrive at reliable error estimates. Tables 5 and 9 contain information on the accuracy with which the chosen LCPs are realized for our simulation parameters. Our estimates for the systematic uncertainties due to deviations from the chosen LCP were then obtained analogously to the case of N f = 2; we refer the reader to appendix B for the details. Here it is worth noting that, similarly to this case, the propagated uncertainties are typically larger than the statistical errors for the R-estimators, eqs. (

Effect of perturbative one-loop improvement
In the lower halves of tables 6 and 7 we give the results for Z A,V after perturbatively subtracting the lattice artifacts to one-loop order. The results have been obtained by first improving the Z A,V determinations for each L/a and g 0 value, and then interpolating to the proper (L 1,2 /a)(β) (see appendix B.5).
Comparing the results for Z A,V before and after perturbative improvement, one sees that the g-definitions are the most affected, and are brought closer to the corresponding l-definitions. All in all, the effect of the perturbative improvement is at most at the level of a couple of percent (cf. figure 7). Hence, not too surprisingly perhaps, the situation is very much the same as for the N f = 2 case.  The individual SF points labelled Z SF A and Z SF A,con are taken from ref. [22] and correspond to the definitions Z A,0 and Z con A,0 , respectively, of that reference. The solid black line is the fit formula to Z SF A also given in [22] and the dashed lines delimit the 1σ region of the fit. Note that this fit function enforces the perturbative 1-loop behaviour for g 2 0 → 0.
In conclusion, our final results for Z A,V are very precise for both LCPs. Similarly to the N f = 2 case, the results for Z A are significantly more accurate than the standard SF determination based on Ward identities [22]. This can be appreciated in figure 6, where the results from table 7 are displayed together with the 2 alternative definitions Z A,0 and Z con A,0 of ref. [22].

Universality and automatic O(a) improvement
Given our estimates for Z A,V we can study the approach to the continuum limit of the ratio between different definitions. We begin with figure 7 where the ratios between the gand l-definitions are considered for the L 1 -and L 2 -LCPs; both the results before and after perturbative improvement are shown. The conclusions are very much the same as for the N f = 2 case. Considering the results before perturbative improvement, along both LCPs, the g and l definitions deviate by at most a couple of per-cent. These differences then perfectly scale with a 2 to zero as the continuum limit is approached. If perturbative improvement is implemented, these differences almost vanish even at the coarsest lattice spacings. There is no significant deviation from a 2 scaling, however, some small admixture of higher order effects cannot be excluded either. It is also interesting to consider the   is also shown. The upper panels show the L 1 -LCP results while the lower ones show those of the L 2 -LCP. In all cases, the dashed lines correspond to linear fits to the data constrained to extrapolate to 1 for a/L 1,2 = 0. Note that the (tiny) effect of the statistical correlation between numerator and denominator has been neglected in these ratios.
continuum limit of the ratio between the definitions belonging to different LCPs i.e. the L 1 -and L 2 -LCP. An example of such a ratio is shown in figure 8. Also in this case, the continuum scaling of this ratio is the one expected, and the initial difference is at the 2 per cent level. Apart from providing an important check of universality and automatic O(a) improvement, these results show that considering one definition or the other for the renormalization of matrix elements of the axial and vector currents, will only introduce small O(a 2 ) differences over the whole range of lattice spacings covered.
a 2 /t 0 X=A X=V Figure 8: Continuum limit of the ratio between the Z l X, L2 definitions, X=A,V, corresponding to the L 2 -LCP, and the Z g X, L1 definitions corresponding to the L 1 -LCP. The dashed lines correspond to linear fits to the data constrained to extrapolate to 1 for a 2 /t 0 = 0.
Finally, we look at ratios between χSF and standard SF determinations. Towards the continuum limit these should also scale like 1 + O(a 2 ), if the SF determinations are O(a) improved. In figure 9 we show the continuum limit of the ratios between the standard SF determinations of ref. [22] and the χSF results of table 7. We here consider both definitions of this reference, and label them as Z SF A = Z A,0 and Z SF A,con = Z con A,0 , respectively (cf. [22] for the exact definitions).
As one can see in figure 9, for their preferred definition, Z SF A , the expected scaling is only setting in around a 2 /t 0 < 0.2, where the SF and χSF determinations differ by a couple of per cent. At the coarsest lattice spacing, corresponding to β = 3.4, the deviation from the O(a 2 ) scaling is significant. The results for Z g A show the largest deviation from the SF determination, which is about 6%. Considering the perturbatively improved χSF results this difference is somewhat reduced to 4-5%, but O(a 2 ) scaling is not observed either. If we consider instead the alternative definition, Z SF A, con , the deviation is reduced to a couple of per cent at the coarsest lattice spacing for Z l A , while, remarkably, the results for Z g A and Z con A,0 are compatible within errors. In particular, the difference between this SF and both our χSF determinations is perfectly compatible with an O(a 2 ) effect over the whole range of lattice spacings considered. While discretization effects can only be defined with respect to some reference definition, we conclude that the alternative SF definition Z SF A,con is, within errors, perfectly scaling with a 2 for β ≥ 3.4 relative to all χSF definitions, whereas the preferred definition Z SF A of ref. [22] requires much finer lattices before this expected asymptotic behaviour sets in. With hindsight, Z SF A,con seems to be a better choice within the SF framework and also has been the preferred SF definition within the N f = 2 setup of ref. [8,42].
Z SF A, con /Z l A, sub Figure 9: Continuum limit of the ratios between the N f = 3 WI determinations of Z A of ref. [22], and the χSF determinations Z g,l A (left panel) and Z g,l A, sub (right panel) of table 7. The Z SF A results are from the fit formula provided in ref. [22], and correspond to their preferred, Z A,0 , definition. The associated dashed lines (red and blue lines) are linear fits to the data with a 2 /t 0 < 0.2, constrained to extrapolate to 1 for a 2 /t 0 = 0. The Z SF A, con results come instead from a fit of the results for the alternative, Z con A,0 , definition considered in ref. [22]. The latter fit was obtained using the same fit ansatz used in ref. [22] for Z A,0 . The associated dashed lines (green and magenta lines) are linear fits to all data, constrained to extrapolate to 1 for a 2 /t 0 = 0.

Summary and conclusions
We have used a new method [20] based on the chirally rotated Schrödinger functional [14] to obtain high precision results for the normalization constants of the Noether currents corresponding to non-singlet chiral and flavour symmetries. The matrix elements of these axial and vector currents play a crucial rôle in various contexts of hadronic physics. Our method differs from the traditional Ward identity method [11,12] in that it compares correlation functions which are related by finite chiral or flavour rotations, rather than infinitesimal ones. The major advantage compared to the Ward identity method consists in the avoidance of 3-and 4-point functions in favour of simple 2-point functions. This very significantly improves on the precision achieved in previous determinations [8,21,22,42,46]. In particular, for the case of Z A , we obtain a reduction of the error by up to an order of magnitude (cf. figure 3 and 6). The relatively poor precision obtained for Z A with the traditional Ward identity methods [8,21,22,42] (around the percent level at the coarsest lattice spacings of interest), has now become a limiting factor in several applications. For this reason, our results are in high demand and have already been used in several works [1,23,47]. In particular, the precise N f = 2+1 scale setting from a linear combination of f K and f π in ref. [1] crucially relies on our values of Z l A in table 6 and the associated uncertainty is negligible compared to the statistical error of the bare hadronic matrix elements. In turn, the precise scale setting result of [1] is entering almost all studies done with CLS gauge configurations: in particular it has enabled the precise result for the 3-flavour QCD Λ-parameter and thus α s (m Z ) by the ALPHA-collaboration [2,39,48,49]. Further applications of our Z A -results include the non-perturbative quark mass renormalization factor in [23] and the related determination of the light and strange quark masses [47]. Regarding the N f = 2 case, the potential improvement of the scale setting in ref. [8] due to our Z A -results would be very significant, too. Tentative estimates anticipate a gain by a factor 3 − 6 in precision, when going from the finest to the coarsest lattice spacing [50].
In order to maximize the usefulness of our results we have chosen the same actions and the same β-values for N f = 2 and N f = 3 lattice QCD as used by the CLS initiative [8,10]. Hence, anyone working with CLS gauge configurations will be able to directly use our results: for N f = 2 we recommend to use Z l A,V, sub from table 4, and for N f = 3 we recommend using Z l A,V, sub either of table 6 or 7. Although the results for Z l A,V, sub are slightly less precise than those for Z g A,V, sub , their L/a-interpolations turn out to be more robust. Furthermore, the effect of the perturbative subtraction of cutoff effects is rather small and only marginally significant with current errors. While the precise choice of the χSF results for the Z-factors is not crucial, it is however very important to be consistent and to not switch definitions when changing β. Only then cutoff effects are guaranteed to vanish smoothly at a rate ∝ a 2 .
Our determination of Z A,V (β) was carried out for each β-value independently, in order to avoid adding statistical correlation between physics results at different lattice spacings. However, if required, it is straightforward to fit our Z-factors to a smooth function of β which interpolates to any intermediate β-value; by matching with perturbation theory at large β one may also obtain results outside the range required for hadronic physics, as done for example in [22]. At some point one then needs to estimate the effect of deviations from the physical box size L = 0.6 − 0.7 fm chosen for our lines of constant physics. However, before this becomes necessary, there is some way to go, as simulations of the χSF up to L/a = 32 (and thus a factor of 2 smaller lattice spacings) are straightforward to perform with current resources.
In applications to hadronic physics some interpolations in β will be necessary if one wants to take into account the O(am) shift in g 2 0 required for consistent O(a) improvement off the chiral limit [26,51]. Fortunately, the corresponding counterterm coefficient, b g , seems to be small, at least in perturbation theory [29], which renders this shift in β negligible at small quark masses.
Looking beyond direct applications of our results in the CLS context, it is quite obvious that the precision gains of this method are generic and could be implemented with any other choice of Wilson type fermions. One would need to implement the χSF boundary conditions following ref. [14], as well as the χSF correlation functions [20]. We also note that the computer resources required are rather modest: in fact our largest lattice size was 16 4 ; indeed, the main work for the present results went into painstakingly following lines of constant physics and the determination of the corresponding uncertainties and their propagation to the Z-factors. We have reported many technical details in the hope that any further applications of the method will be able to benefit from our experience. One possible improvement we did not explore was to measure the derivatives (3.9,3.10) by computing the corresponding operator insertions into the correlation functions directly on the tuned ensembles; this was e.g. done in refs. [1,2] for the PCAC mass, t 0 , and other observables, and this would certainly allow one to further improve on the precision, as no assumptions on the derivatives need to be made.
Future applications of methods based on the χSF include the determination of the ratio between scalar and pseudoscalar renormalization constants, Z S /Z P . Advantages of the χSF are also expected for scale-dependent problems, such as the renormalization of 4quark operators, where the contamination by O(a) effects could be significantly reduced by the mechanism of automatic O(a) improvement [19]. Finally the χSF offers new methods for the determination of O(a) improvement coefficients, which we hope to explore in the future.

Acknowledgments
We would like to thank Rainer Sommer and Stefano Lottini for useful discussions, and Patrick Fritzsch and Tim Harris for helpful comments. We are moreover thankful to the computer centers at ICHEC, LRZ (project id pr84mi), and DESY-Zeuthen for the allocated computer resources and support. The simulation code we used for the simulations is based on the openQCD package developed at CERN [36].

A Simulation parameters and results
L/a β κ z f mL × 10 3 g ud A (L/2) × 10 3 P Q N ms  Parameters of the N f = 2 ensembles and corresponding results for mL and g ud A (L/2). The total number of measurements we collected is given by N ms ; these are spaced by 10 MDUs. In the table we also give the percentage P Q of gauge fields with Q = 0.  Parameters of the N f = 3 ensembles and corresponding results for mL and g ud A (L/2). The total number of measurements we collected is given by N ms ; these are spaced by 10 MDUs. In the table we also give the percentage P Q of gauge fields with Q = 0.    table 10: we thus obtain the same values. The systematic errors associated with maintaining the condition (4.2), instead, involve different derivatives for the data in table 10 and the perturbatively improved ones (cf. appendix B).  Renormalization constants Z g,l A,V for the different N

B Error propagation, systematic error estimates and comparison with perturbation theory
In this appendix we describe in some detail the elements required to carry out the steps 2,3 and 4 sketched in subsection 3.6, for the propagation of uncertainties. We first report on the numerical estimates of the various derivatives (steps 2,3) for both N f = 2 and N f = 3. These estimates are obtained on smaller lattices L/a = 8 and L/a = 6, 8 respectively, and then used for all lattice sizes. We therefore also summarize the expected scaling with L/a of these derivatives and confirm this to first non-trivial order in perturbation theory. Finally, the interpolation of the Z-factors in L/a (step 4, where necessary) is discussed, first for N f = 2, where this step is almost avoidable, and then for N f = 3, where interpolations are required at most β-values, thus requiring a more thorough analysis.
B.1 Estimating the derivatives: N f = 2 As described in section 3.6, to estimate the uncertainties in am cr and z * f we require the derivatives (3.9) and (3.10). We estimated these through dedicated simulations at L/a = 8 and β = 5. We observe that the z f -derivative of mL vanishes within an uncertainty much smaller than the values of the other derivatives, so that we may safely set it to zero. These results were then used for all other ensembles and lattice sizes listed in table 8. The uncertainty in the tuning of mL and g ud A (L/2) is thus converted to one in m cr L and z * f . Specifically, one has that, and ∆(mL), ∆g ud A , ∆(m cr L) and ∆z * f , are the uncertainties in mL, g ud A (L/2), m cr L and z * f , respectively. For the uncertainties ∆(mL) and ∆g ud A we took the (absolute) values of mL and g ud A (L/2) measured on the given ensemble (cf. table 8), plus 2 times the corresponding statistical errors. The corresponding systematic errors on the Z-factors can then be estimated as, Through the very same simulations used to determine the derivatives (B.1) we obtained, For N f = 3 we proceeded in very much the same way, except that we carried out a more complete study of L/a = 8 lattices at all β-values, except β = 3.85, and also included additional L/a = 6 lattices at β = 3.4 and 3.46. As before for the derivatives (3.9) we used their mean values and set dmL/dz f = 0, while for the derivatives (3.10) we considered their (absolute) mean values plus twice their statistical errors. However, here we did this separately for all β-values, with β = 3.7 results also applied at β = 3.85. Table 12 contains the results for Z g,l A,V for all ensembles of table 9, including the different estimated uncertainties. As with N f = 2, the soundness of our assumption is supported by experience gained during the tuning runs to find m cr and z * f .

B.3.1 Expected L/a-scaling
Using general arguments based on the Symanzik expansion and P 5 parity [20] one may obtain the expected scaling with the lattice spacing a of the derivatives of the PCAC mass and g ud A with respect to m 0 and z f , and To explain how one arrives at these scaling properties we go through the example of the PCAC mass: where we have introduced the modified PCAC mass m through the equatioñ ∂ 0 g ud A;z f = 2m g ud P;z f , (B.9) and the notation ; z f indicates differentiation with respect to z f [20]. The crucial point to note is that this differentiation merely modifies the fields at the time boundaries and thus produces a different matrix element for the PCAC relation and hence the modified PCAC mass m . Since the difference m − m between two PCAC masses is of O(a) in general, and the z f -derivative of the P 5 -even correlation function g ud P is P 5 -odd and thus of O(a), we arrive at O(a 2 ) for the complete expression. It is re-assuring to see that this derivative is indeed found to be small in the simulations.
Similar arguments lead to where the derivatives are taken at fixed β, z f and β, am 0 , respectively.

B.3.2 Comparison with perturbation theory
Perturbation theory 5 confirms all of these expected scaling properties. The derivatives (3.9) have only been considered to tree-level, which gives, , For the mass sensitivity of the axial current we have, to one-loop order, and, for the vector current, Here, the one-loop coefficients are the values at L/a = 12 and are stable within 3-10 per cent for the range L/a from 8 to 16. Incidentally, this result resolves qualitatively a puzzle posed by the non-perturbative results, eq. (B.4), where the derivatives of the gand ldefinitions are almost the same in magnitude and opposite in sign, whereas the tree level results are 1 and 0, respectively, as first noticed in [20].
The z f -sensitivity of the Z-factors is easily described: all z f -derivatives vanish at tree level and the one-loop coefficients for the g-defintions are very small and vanish with a rate roughly proportional to a 3 for Z g A,V and both gauge actions. On the other hand, the l-definitions behave as expected (s. above): very similar numbers are obtained which are, for both gauge actions, within a few percent given by For completeness we report the perturbative results for am cr and z * f we obtain from the same computation. For the critical mass to order g 2 0 , the known values [33,44,52,53] (with C F = 4/3 for gauge group SU(3)), are already reproduced to 4-5 digits on lattices with L/a in the range from 8 to 16. As for z * f , we have, to order g 2 0 and for a/L → 0, (1) where the Wilson action result is from ref. [20], whereas the LW action value is the one of L/a = 16 with a generous guess for the error. Also in this case the values at finite L/a from 8 and 16 coincide with these numbers to 4-5 digits precision. We observe that the quantitative comparison of non-perturbative data with bare perturbation theory to order g 2 0 works quite well in certain cases. For instance, for N f = 2 at β = 5.2, we compare the non-perturbative value am cr = −0.3282 to am cr = −0.3116 + O(g 4 0 ), and similarly for z * f = 1.288 we need to compare to z * f = 1.258 + O(g 4 0 ). Also the non-perturbative current normalization constants themselves are reproduced by one-loop perturbation theory at the 5-10 percent level (compare table 1 with tables 4 and 6,  7). On the other hand, the majority of the derivatives differ very significantly, for instance, is off by a factor 2, and for the l-definition the comparison is between −0.127(9) and −0.023, which differs by a factor 5. For the z f -derivatives we note that perturbation theory correctly predicts the smallness of the sensitivity in the g-definition. Quantitatively, the perturbative z f -derivatives of Z l A,V at β = 5.2 and L/a = 8 are −0.046 + O(g 4 0 ), to be compared with eq. (B.5), so again we observe a difference by a factor 2.
Finally, for the interpolations in L/a of the data (s. below) we are also interested in the derivatives of the Z-factors with respect to x = a 2 /L 2 at fixed bare coupling. In perturbation theory we can obtain approximate results from table 1. Expanding the x-derivatives to order g 2 0 are approximately given by provided higher order cutoff effects are small. We find that this is quite well satisfied, with very similar coefficients for both Wilson and LW actions, given approximately by whereas k l A,V are ca. 20 to 30 times smaller in magnitude, for axial and vector cases, respectively, and come with the opposite sign. Comparing this with the β = 5.2 nonperturbative data in Eqs. (B.23) we see again that for the g-definitions these derivatives are reproduced by perturbation theory up to a factor 2, while for the l-definitions, perturbation theory to O(g 2 0 ) is clearly missing the bulk of the effect. While, as expected, the nonperturbative derivatives are smaller than for the g-definitions, it seems that the smallness of the O(g 2 0 ) term is an accident and higher orders are dominating at these values of β. To conclude this comparison, perturbation theory often gives valuable qualitative information and may provide reasonable starting values for the tuning of am 0 and z f . However, quantitatively, the agreement with non-perturbative data at lattice spacings of interest for hadronic physics hugely varies for different observables. Hence, the main practical use of perturbation theory consists in the perturbative subtraction of cutoff effects. Here, even a qualitative agreement, which may be quantitatively off by a factor 2, still means a welcome reduction of cutoff effects by 50 percent, and our data analysis does indeed point to such benefits.
Given this situation, we have refrained from using perturbative data in our estimates of the derivatives, and we have decided to ignore the favourable O(a) scaling, eq. (B.10), when applying the results obtained at L/a = 8 at all other L/a-values, too. In this respect, the tuning runs for am 0 and z f , both for N f = 2 and N f = 3, provided some consistency checks which make us confident that the chosen procedure is indeed sound and rather conservative.
B.4 Interpolation in L/a for N f = 2 Once the uncertainties associated with the conditions (2.13) have been propagated to Z A,V at given β and L/a, one needs to keep to the LCP condition (4.2) and also take into account the corresponding uncertainties. The choice (4.2) is made such that the L/a = 8 results at β = 5.2 satisfy this condition by definition, while at β = 5.3 we needed to interpolate the results for L/a = 8, 10, 12 to the target value (L/a)(5.3) = 9.18(21) (cf. table 3). We have performed a simple linear interpolation in (a/L) 2 which describes the data very well, similarly to the case N f = 3 which will be discussed in more detail below. For β = 5.5 and 5.7, the simulated L/a values are, within errors, compatible with the target values in table 3. Systematic errors related to the condition (4.2) were then estimated by using the slope of the interpolation in x = (a/L) 2 at β = 5.3 also for the higher β-values. Note that we interpolate results with bare parameters am cr and z * f tuned at the given βand L/a-values. Hence the derivative defines the sensitivity to a change of the physical size of the system. This is a pure cutoff effect of O(a 2 ) on the Z-factors. Since, by the choice of the variable x, a factor a 2 is also divided out, the x-derivative is expected to be of O(1) and we expect a smooth dependence of this derivative on β; this expectation is in fact confirmed by the results for N f = 3 where the slope shows a very mild β-dependence over the whole range. As we are looking at an O(a 2 ) effect in disguise, it is no surprise that the results depend on whether or not the cutoff effects have been subtracted perturbatively.
Without perturbative subtraction we obtained the results, The corresponding systematic error is then simply taken to be, where x(β) = ((a/L)(β)) 2 and σ(x(β)) is the associated statistical error. As a further safeguard we took for the derivatives (B.23) and (B.24) the (absolute) mean value plus twice their statistical error. We observe that the l-definitions have a milder L/a-dependence than the g-based ones, unless perturbative improvement is implemented.
B.5 Interpolation in L/a for N f = 3 Once the systematic errors deriving from the tuning of am 0 and z f have been taken into account, the results for Z A,V at different L/a and fixed β must be interpolated to either (L 1 /a)(β) or (L 2 /a)(β), depending on the LCP; for completeness the values of Z A,V prior interpolation are given in table 12. We have considered three types of interpolation in x = (a/L) 2 , these are: linear using all 4 available values of L/a (cf. table 5), linear using only the 3 closest L/a-values to the target L 1,2 (β), and quadratic using all 4 L/a-values. Given this choice, the interpolations needed for the L 1 -and L 2 -LCPs only differ for β = 3.4, where, by definition, L 1 /a = 8 is exact, and in the case of linear interpolations with 3 points at β = 3.55. Recall that for β = 3.85 no interpolation is required, as L 2 /a = 16 is exact and this β-value has been excluded for the L 1 -LCP. Starting with the L 1 -LCP, the different interpolations describe the data quite well in general, particularly so for the results at the two smallest lattice spacings and for definitions based on the l-correlators. The most relevant exception is given indeed by the linear interpolation of Z g V at β = 3.46 using all 4 values of L/a, for which we find a χ 2 /d.o.f ≈ 2.2. It should be noted, however, that the χ 2 -criterion does not come with the usual probability interpretation due to the errors being dominated by systematics. In any case, the interpolated values are generally compatible at the 1σ level. Considering the perturbatively improved data, the quality of the interpolations is generally improved, and all fits have excellent χ 2 . The beneficial effect of the perturbative improvement can be appreciated by comparing figure 10 and 11, where the Z A,V interpolations at β = 3.46 are shown for the cases before and after perturbative improvement, respectively. This example also illustrates the general feature that, before perturbative improvement, the l-definitions have a significantly milder L/a-dependence. Based on these observations, we take as our final estimates for the Z-factors the results of the quadratic fits, which have the largest errors.
Regarding the L 2 -LCP, the situation is more complicated due to the fact that we need to interpolate the data at the coarsest lattice spacing, β = 3.4. The quality of the interpolations is still good in general, but there are a few significant exceptions. We note that all these cases involve g-definitions: indeed, we obtain a pretty large χ 2 /d.o.f. for the linear fits of Z g V at β = 3.4 and 3.46, around 7.8 and 5 respectively. Also the quadratic fit for Z g A at β = 3.4 has a large χ 2 /d.o.f ≈ 1.8. While in this case, however, the results of the interpolation are compatible with those of the linear fits within less than one standard deviation, in the case of Z g V at β = 3.4 the discrepancy between the linear and quadratic interpolations is close to 3 standard deviations. The situation definitely improves when the perturbatively improved data are considered. In this case, with the exception of the linear interpolations with 4 points of Z g V,A;sub at β = 3.4, all fits have very good χ 2 , and give compatible results within one standard deviation or so. As in the case of the L 1 -LCP, we take as our final estimates for Z A,V the results of the quadratic fits, which have the best χ 2 -values and the largest errors.