Non-perturbative determination of c_V, Z_V and Z_S/Z_P in N_f=3 lattice QCD

We report on non-perturbative computations of the improvement coefficient c_V and the renormalization factor Z_V of the vector current in three-flavour O(a) improved lattice QCD with Wilson quarks and tree-level Symanzik improved gauge action. To reduce finite quark mass effects, our improvement and normalization conditions exploit massive chiral Ward identities formulated in the Schroedinger functional setup, which also allow deriving a new method to extract the ratio Z_S/Z_P of scalar to pseudoscalar renormalization constants. We present preliminary results of a numerical evaluation of Z_V and c_V along a line of constant physics with gauge couplings corresponding to lattice spacings of about 0.09 fm and below, relevant for phenomenological applications.


Introduction
A popular discretization for quark fields on the lattice are Wilson fermions. However, as a consequence of removing the unwanted doublers in the naive lattice fermion action, it exhibits leading cutoff effects of O(a) and the explicit breaking of chiral symmetry. As for the former, a systematic way of resolving this is the Symanzik improvement programme, which amounts to add the so-called clover term to the action and further irrelevant operators to local composite fields, canceling their O(a) corrections, while the latter is accounted for by introducing finite renormalization constants. To eliminate all O(a) contributions from physical quantities and to restore chiral symmetry at this order, these improvement counterterms and renormalization factors have to be fixed non-perturbatively.
In this work we specifically look at the renormalized and improved isovector current, which in the chiral limit of vanishing sea quark masses and at non-zero valence quark mass can be parametrized as with (V I ) a µ (x) = V a µ (x) + ac V∂ν T a µν (x) (2) = ψ(x)γ µ τ a 2 ψ(x) + iac V∂ν ψ(x)σ µν τ a 2 ψ(x) , where τ a acts in flavour space and∂ µ denotes the symmetric lattice derivative. The quark mass dependent O(a) improvement term proportional to b V , which corrects for quark mass dependent cutoff effects, was recently calculated non-perturbatively for N f = 3 in [1,2]. The renormalization constant Z V and the improvement factor c V , however, have not yet been investigated non-perturbatively so far in the case of three-flavour QCD and are subject to this work. Potential applications of the vector current and its matrix elements include computations of semi-leptonic decay form factors and of the timelike pion form factor, as well as contributions to the anomalous magnetic moment of the muon and thermal correlators related to the di-lepton production rate in the quark-gluon plasma. Another object of interest is the ratio Z S /Z P of the scalar and the pseudoscalar renormalization constants, which plays a rôle in relating renormalized PCAC and subtracted quark masses to each other. The two constants themselves exhibit a scale dependence that cancels in the ratio, though. Accordingly, in the zero sea quark mass limit and at non-vanishing valence quark mass, the corresponding renormalized currents are defined as and already comply with O(a) improvement without any correction terms.

Renormalization and improvement conditions
All improvement and renormalization conditions explained below involve the O(a) improved PCAC quark mass defined as with standard notation for (symmetric, backward and forward) lattice derivatives and where c A for O(a) improved N f = 3 lattice QCD with Wilson fermions [3] and tree-level improved gauge action, as employed here, is non-perturbatively known from [4]. Let us mention that there exists a very promising alternative approach to determine renormalization factors through imposing appropriate conditions based on the PCAC relation in the Schrödinger functional with chirally rotated boundary conditions, see [5], where it was also tested in perturbation theory. Apart from its advantage of entailing automatic O(a) improvement, it turned out that, e.g., in case of the renormalization factor of the axial current for N f = 2, more precise results than with standard Schrödinger functional boundary conditions can be obtained [6].

Renormalization of the vector current
The renormalization condition for the vector current is derived from the vector Ward identity [7], By choosing the spacetime region R to consist of all times smaller than x 0 , the only contribution stems from the timeslice x 0 which results in We identify the operators O int and O ext with the boundary fields at x 0 = 0 and x 0 = T , where ζ and ζ are the Schrödinger functional boundary fields at x 0 = 0 and their primed versions the fields at x 0 = T , respectively. After replacing both sides by their renormalized lattice counterparts, we arrive at with The Ward identity is valid for all x 0 , although boundary effects are expected far from the temporal center of the lattice. In order to get a better handle on statistical fluctuations, we have evaluated the renormalization condition at the central four timeslices and taken the average.

Improvement of the vector current
The improvement condition for the vector current was first presented in [8] and is based on the axial Ward identity. By insertion of an axial current as an operator inside the spacetime region R we get By specifying R as the region between the timeslices x 0 = t 1 and x 0 = t 2 with t 1 < y 0 < t 2 , two surface terms arise: Since the Ward identity is valid for all ν, we take ν = k. The source operator O c ext is chosen as where ζ and ζ are quark fields at the boundary x 0 = 0. After implementing this improvement condition in terms of Schrödinger functional correlation functions, one finds , (14) omitting the sea and valence quark mass b-coefficients for brevity; it is to be understood as our final expression, which can be solved for c V (once Z V and Z A are known). For explicit definitions of the correlators we refer, e.g., to [7][8][9][10], with contributions that are diagrammatically represented via quark diagrams corresponding to possible Wick contractions in figure 1; more details will be given elsewhere [11]. For the present analysis, (14) was evaluated at t 1 = T/4 and t 2 = 3T/4, as originally suggested in [8]. Considered as a function of the timeslice variable y 0 , a plateau at the temporal center of the lattice is identified for the (local) c V (y 0 ). In order to tame statistical fluctuations, the quoted preliminary values for c V are extracted as averages of the central two timeslices.

Ratio of renormalization constants Z S /Z P
To derive a renormalization condition for the ratio of quark mass renormalization factors Z S and Z P , we exploit a renormalization condition that once more is derived from the massive axial Ward identity -closely following the ALPHA Collaboration's method to compute Z A for N f = 2, 3 [9, 10] -, but now relying on a pseudoscalar insertion as internal operator with a behaviour under variation: For d abc not to vanish (and thus to ensure sensitivity to the scalar density on the r.h.s. of this equation), one has to work with a SU(N f ) algebra in the valence sector, where N f ≥ 3. Here, we choose SU (3), assume a b and adopt a product of two pseudoscalar boundary sources, compare (8), as external operator in the integrated axial Ward identity which is similar to (11) but involves a pseudoscalar density insertion with a variation according to (15). Upon identifying each piece with a Schrödinger functional correlator, some steps of algebra [12] yield a formula that can be solved for Z P /Z S (once Z A is known) and in which the intrinsic scale dependence of the individual renormalization factors drops out, viz.
Again, any b-coefficients are suppressed here. The correlators are defined analogously to those in [7][8][9][10], and their explicit forms will be given fully elsewhere [12].   (18). Filled (open) circles stand for the creation (annihilation) of a quark at the boundaries of the lattice, while squares indicate insertions of local composite fields. In a first trial analysis, we had evaluated (18) for insertion points y 0 = T/2, t 1 = T/3 and t 2 = 2T/3. Still, we do not quote any results for Z P /Z S in this status report, because a careful study to extract it from the renormalization condition proposed here has only started after the conference.

Simulation details
As improvement coefficients and renormalization factors are short-distance quantities, they can be extracted by imposing suitable conditions in a finite (i.e., in practice, small) physical volume. This is realized by the Schrödinger functional framework, governed by periodic boundary conditions in space and Dirichlet ones in time. The gauge field configuration ensembles used in this work are almost identical to the ones that were generated in the context of the improvement and renormalization of the axial vector current [4,10] and cover the β-range of the N f = 3 large-volume QCD configurations of the CLS effort, corresponding to lattice spacings of about (0.05 a 0.09) fm [13,14]. Their specifications are collected in table 1. To supplement the data base of configurations already available from [4,10], the production of a few new ensembles -labeled by A1k3, D1k2 and D1k3 in the table -was started. These ensembles exhibit a more chiral (i.e., closer to zero) mass of the three mass-degenerate sea quarks and thereby allow for getting a better handle on the mass dependencies of the quantities of interest that will prove to be essential in the case of c V . Compared to the previous N f = 0 study [8], we have implemented various refinements: First of all, as detailed in [4], all gauge field ensembles entering the analysis lie on a line of constant physics characterized by a fixed spatial physical volume of L ≈ 1.2fm = constant, T ≈ 3L/2 and almost vanishing mass of the (degenerate) sea quarks. The valence quark mass in the computation of correlation functions equals the sea quark value. This entails that the renormalization and improvement factors become smooth functions of the bare coupling, i.e., g 2 0 = 6/β. Only the ensemble B2k1 deliberately misses the condition of fixed physical volume and is used to quantitatively investigate the effect of such a deviation on the results. Furthermore, again following [4], the Schrödinger functional correlation functions incorporate optimized boundary wave-functions, in order to suppress excited state effects and thus to maximize the overlap with the ground state in their spectral decomposition. Finally, for the case of c V , we also have identified the importance of the additional mass term in the axial Ward identity, (14), the impact of which will be discussed in the results section. In this context, we have tested different sets of insertion times for the individual operators and found a specific choice that seems to reduce the effects caused by the non-zero mass comprehensively.
The statistical error analysis of the Markov chain Monte Carlo data utilizes the Γ-method based on evaluating autocorrelation functions [15] and was cross-checked against binned Jackknife estimates.

Results
The analysis underlying the results presented here was done including all topological sectors. By virtue of the theoretical argument that our results -being based on Ward identities, as operator identities holding in any topological sector -should be insensitive to the topological charge Q top , we believe that the influence of restricting the computations to one sector of fixed Q top (say, Q top = 0) is negligible, modulo the accompanying reduction of statistics. This expectation still needs to be confirmed in the final analysis, though.
The left panel of figure 2 shows an representative evaluation of the local PCAC quark mass, am PCAC (x 0 ), on the gauge configurations of ensemble C1k3. The actual values for all PCAC masses entering our analysis are always chosen as plateau averages over the central L/2 timeslices of the temporal extent of the lattice. Thereby it is guaranteed that also these averaging intervals are scaled in physical units in the same way as all other length scales, in order to obey the constant physics condition in all steps of the computation.

Z V
Results for the vector renormalization constant Z V are presented in the right panel of figure 2, in comparison to 1-loop perturbation theory taken from [16]. The individual data points, obtained by a chiral extrapolation to zero (valence = sea) quark mass using b V in (9) from [1, 2] (but neglecting the corresponding b-coefficient in the sea quark sector), show a smooth behavior that is well described by a polynomial fit constrained by perturbation theory. The preliminary interpolation formula reads:

c V
As outlined above, the condition of [8] to fix the vector current improvement coefficient c V is extended by accounting for an additional term that naturally arises when the Ward identity is employed at finite quark mass, cf. (14). The impact of this term is demonstrated in the left panel of figure 3, where the chiral extrapolations (using values for the associated valence quark mass b-coefficients from [1,2]) for L/a = 16 and g 2 0 = 1.7084 (β = 3.512) with and without the mass term are compared with each other. Although both extrapolations nicely meet in almost the same chiral limit at am PCAC = 0, the data without inclusion of the mass term show a much steeper behaviour. This finally results in a larger error at am PCAC = 0 and thus underlines the importance of refining the improvement condition for c V through accounting for the mass term in the analysis even at small but finite quark masses.  [17]. The gray data point refers to a tentative analysis on the gauge configuration ensembles D1k2 and D1k3, whose generation was launched during the conference. Hence, it is only indicative and not yet included in the fit.
Additionally, to quantitatively check for the influence of a violation of the constant physics condition on our analysis, results from the gauge configuration ensemble B2k1 are displayed in the same figure. This ensemble (see table 1) has a β-value shifted by an amount, which corresponds to a ∼ 6% shift in the spatial extent of the physical volume, and thus induces a significant deviation from our condition L ≈ 1.2fm = constant. As can be seen in the left panel of figure 3, these data points align fairly well with the other points along the (linear) fit function. Hence, we conclude that any deviations from the constant physics condition of this order of magnitude or below are of only minor influence and can safely be neglected on the level of the final precision for c V . 1 Note that this is also in line with the findings already reported in [4,10].
The preliminary estimates for c V are presented in the right panel of figure 3, together with the prediction from 1-loop perturbation theory that we have extracted for our lattice action from the perturbative results in [17]. The gray data point stems from the ensembles D1k2 and D1k3 (see table 1). Since the generation of these gauge field configurations was only started during the conference and is still ongoing, we exclude it from the subsequent analysis steps for the purpose of the present status report. Nevertheless it is reassuring that this -so far only indicative -result appears to blend in well with the g 2 0 -dependence of the other points. At this point, we therefore describe our results for c V by a preliminary interpolating Padé fit, constrained by 1-loop perturbation theory in the asymptotic g 2 0 → 0 regime, as where, as stressed above, the gray point in the right plot of figure 3 is not included in the fit. Moreover, any uncertainties originating from Z V or Z A (entering the final formula for c V according to (14)) have not yet been propagated into the errors on c V quoted in the figure such that we still expect them to slightly increase after a final analysis. Note that, in qualitative agreement with observations already made in the exploratory quenched study [8], the non-perturbative c V substantially deviates from perturbation theory in the range of bare couplings (resp. β-values) typically encountered in large-volume applications with the lattice action employed here.

Outlook
For the completion of our work to determine the renormalization and improvement factors discussed in this report it essentially remains to 1.) evaluate the relevant correlators for the full statistics on all ensembles of table 1, 2.) check for independence of the results on topology by repeating the computations in the sectors of fixed Q top , 3.) quantify the size of possible O(a) ambiguities in improvement (resp. renormalization) conditions for the vector current and, in particular, 4.) to also perform the data analysis to extract the ratio Z S /Z P . For a related study to calculate improvement b-coefficients in the valence sector, multiplying mass dependent O(a) Symanzik counterterms to local operators, as well as the ratio Z m Z P /Z A of quark mass renormalization constants, see [18].