Mass-improvement of the vector current in three-flavor QCD

We determine two improvement coefficients which are relevant to cancel mass-dependent cutoff effects in correlation functions with operator insertions of the non-singlet local QCD vector current. This determination is based on degenerate three-flavor QCD simulations of non-perturbatively O(a) improved Wilson fermions with tree-level improved gauge action. Employing a very robust strategy that has been pioneered in the quenched approximation leads to an accurate estimate of a counterterm cancelling dynamical quark cutoff effects linear in the trace of the quark mass matrix. To our knowledge this is the first time that such an effect has been determined systematically with large significance.


Prelude
Lattice QCD has reached a maturity in which some quantities can be computed very precisely such that QED or isospin breaking effects become relevant when compared to experiment. Other quantities are being improved in various ways to match the accuracy reached in state-of-the art experiments, or to use it as input in the search for new physics.
In either case, the control over continuum extrapolations of lattice data is an important ingredient that contributes to the overall precision in any observable. If the QCD action is discretized using Wilson fermions [1], chiral symmetry is broken explicitly by a term proportional to the lattice spacing a, and consequently (renormalized) observables receive contributions linear in the lattice spacing. With some additional effort these terms can be systematically removed as outlined by Symanzik [2][3][4]. In the present paper we aim at a computation of the coefficients b V andb V which multiply mass-dependent counterterms to restore a continuum scaling to O(a 2 ) in the local vector current made of massive quarks. Thus they are relevant in calculations such as semileptonic vector form factors or the muon anomalous magnetic moment.
In cases where the impact of an improvement coefficient seems negligible from a perturbative point of view, it can still be important to have a non-perturbative estimate at hand in order to properly judge at what level of accuracy they will become important.
To determine such coefficients one requires well-justified improvement conditions. These conditions can be chosen freely to some extent but differ in complexity and quality regarding statistical and/or systematic effects-like ambiguities at higher order in a. They all remove the leading scaling violation they have been designed for, sometimes at the cost of introducing large contributions at higher orders. The condition explored in the present paper relies on the QCD isospin vector symmetry, and we argue that this method, adopted from reference [5], has little systematic effects. Being implemented in the QCD Schrödinger functional (SF) [6][7][8][9] it furthermore allows to efficiently explore the region about the chiral limit over a wide range of couplings and masses at relatively small computational costs. As a result we are able to differentiate between the valence and dynamical (sea) quark sector when the gauge coupling becomes strong.
The pattern of chiral symmetry breaking, inherent to Wilson fermions at finite lattice spacing, intertwines renormalization and O(a) improvement [10]. As a result, the local (isovector) vector current, V a µ (x) = ψ(x)γ µ 1 2 τ a ψ(x), of mass-degenerate fermions renormalizes as 1 in a massless renormalization scheme [5]. Here, V I is the bare vector current, on-shell O(a) improved through the symmetric lattice derivative of the tensor current T a µν with appropriate improvement coefficient c V . We follow the standard notation in the literature [10,11] where c * terms cancel O(a) cutoff effects and b * ,b * the leading mass-dependent effects of order am, induced from the valence or sea quark sector, respectively. 2 These coefficients are considered functions of the bare gauge coupling g 2 0 while the finite renormalization Z V has to be evaluated atg 2 [10,13,14] to maintain full O(a) improvement of V R in the presence of massive quarks. Rather being a mass-dependent finite renormalization than a genuine lattice artefact, the coupling counterterm coefficient b g (g 2 0 ) is required to cancel the mass-dependence of g 2 0 in such a way thatg 2 0 is mass-independent up to terms of order a 2 . Note that both couplings coincide at the critical line. From eq. (1.1) it then becomes clear that b V andb V can be determined by individual varying the valence quark mass am q and sea quark mass Tr [aM q ].

Centerpiece
We are simulating N f = 3 mass-degenerate flavors of Wilson fermions with tree-level improved gauge action in the Schrödinger functional, exhibiting Dirichlet boundary conditions in time and periodic boundary conditions in spatial directions. For on-shell O(a) improvement of the action we use the non-perturbatively determined clover coefficient c sw [15] of ref. [16] and one-loop values of the boundary counterterms c t andc t [10,17,18]. Except for simulating at non-vanishing quark mass, our setup thus equals the one of ref. [18] with vanishing boundary fields, T = L and θ = 1/2. In the following we adopt the notation of refs. [5,10,19].

Improvement condition of the vector current
For completeness we first introduce the matrix elements used in the present calculation and comment on possible systematic effects. In the SF sources are typically defined on the time boundaries at Euclidean times x 0 = 0 and x 0 = T , respectively. They are projected to zero momentum and transform according to the vector representation of the exact isospin symmetry. The simplest correlation function is the two-point boundary-to-boundary correlator which up to terms of order a 2 equals the three-point boundary-to-boundary correlator with an insertion of the time component of the renormalized vector current, eq. (1.1).
It is important to notice that (a) the normalization and improvement of the sources are irrelevant and (b) the O(a) counterterm proportional to c V cancels as well. Thus two potential sources to O(a 2 ) ambiguities are absent and the condition becomes and Z V the overall normalization in the chiral limit. It also implies that (c) the relation is independent of x 0 up to cutoff effects, or more specifically, independent of c t andc t . However, if not stated otherwise we take x 0 = T /2 as the preferred definition of relation (2.4). Our restricted interest in the mass-improvement coefficients leads to another simplification. Estimators of b V andb V can be realised as derivatives such that (d) the overall normalization factor Z V cancels, e.g., The derivative essentially eliminates all mass-independent cutoff effects such that R V = b V +O(am) mainly carries ambiguities that depend on the quark mass. 3 If the derivative is computed in a range of quark masses where f 1 /f V is dominantly linear, these ambiguities are practically absent as (e) it does not matter at which point the derivative is exactly evaluated. As a consequence of (d)+(e) we conclude that one does not necessarily need to impose a line of constant physics when determining R V . Concerning sea quark mass effects we remark that in the unitary case (m q = m sea ) with N f degenerate quarks ([M q ] ij = δ ij m sea with i, j = 1, . . . , N f ) only the combination b V + N fbV is accessible from eq. (2.4). We thus definē at M q = 0 and g 2 0 fixed, (2.6) as estimator for b V + N fbV and obtainb V trivially from an appropriate linear combination with R V . This being said, we proceed as follows: (1) at different values of the lattice spacing we use a given set of dynamical three-flavor simulations with varying sea quark mass to determine the chiral limit (m sea = 0) by simple interpolation, (2) we determine the combination b V + N fbV for the unitary case via eq. (2.6), (3) on the ensemble in the vicinity of vanishing quark mass we vary the valence quark masses in order to determine the mass-derivative (2.5) and thus b V . With above arguments about the potential quality of the improvement conditions (2.5,2.6), and to keep the numerical effort and costs affordable, we restrict our simulations to lattices of size L/a = 8. We show explicitly that this is sufficient in the present case. Our data analysis accounts for autocorrelations via the Γ-method [20] and correlations are included through standard linear error propagation.

Chiral trajectory and current quark mass in finite volume
As usual we define the chiral limit as the point where the current quark mass m vanishes. 4 For simplicity we restrict the discussion to the unitary case (m q = m sea ) and remark that overall (multiplicative) renormalization factors are not required. In the O(a) improved theory, the relation between current quark mass m and bare subtracted sea quark mass m sea reads whereZ = Zr m is a relative (finite) normalization and B a well-defined combination of mass-improvement coefficients [11]. With similar numerical simulations close to the chiral limit at hand [23], we do apply one iteration of refinements at 9 different values of the gauge coupling β = 6/g 2 0 to perform simulations in close proximity of m = 0. In all cases we are as close as |am| < 0.001 and aim for a symmetric variation in the quark mass spanning several orders of magnitude to probe the regime of linearity of eq. (2.7) with significant precision. Where it is possible we go as far as Lm ≈ ±2, but only compile the subset most relevant for the determination of is the bare subtracted quark mass based on the critical hopping parameters κ crit of table 1. By steadily increasing the number of points about vanishing current quark mass entering a least square minimisation to a straight line in 1/(2κ sea ), we inspect the quality of the fit to decide on the range of linearity. In all cases, we do not find relevant deviations in the range |am sea | < 0.014. As this is of minor importance anyway, we determine κ crit by fitting to data in that range only. For completeness we give the corresponding slope,Z, as estimate for the potential size of Zr m , but have to remind the reader that this quantity is better determined along a line of constant physics. We present a typical interpolation in the top panel of figure 1, compiled from data of table 2 at an intermediate coupling (β = 4.0) and strong coupling (β = 3.4). For the latter we know the size of the lattice spacing, a ≈ 0.086 fm, from ref. [24]. The bare subtracted quark mass at am sea = 0.1 thus corresponds to m sea ≈ 230 MeV, and at am sea = 0.01 accordingly to 23 MeV. To state the obvious: within the given uncertainty it is confirmed that the improvement B-term in eq. (2.7) is negligible for up-/down-like quarks while it can be relevant for strange-like quarks, and certainly is for heavier ones.  The inner graph is a magnification of data entering the straight line fit close to the chiral limit (|am q | < 0.014), while the outer graph visualises the deviation from linear behaviour far away from the chiral limit. All data points are independent from each other and taken from table 2. The quality of fit tends to be somewhat better for Z V data when compared to Lm data. Note that the fit quality of Lm at β = 3.4 improves when the points atop |am q | = 0.012 are included.

Unitary quark mass dependence
We have measured the correlation functions f 1 and f V in the unitary setup for all sea quark masses. Their values are listed in table 2 in dependence of m sea and we use Z V (m q ) as short hand for the mass-dependent ratio f 1 (m q )/f V (m q ). We clearly profit from correlations between the two boundary-to-boundary correlation functions but note that values on different ensembles are uncorrelated. As in the previous section we study the range in which the data is well compatible with a straight line fit ansatz, which again is the case for |am sea | < 0.014, cf. bottom panel of figure 1. At each coupling β we neglect the ensemble closest to the chiral limit in order to be fully independent from the determination of R V presented in the next section. The associated interpolation provides the intercept at m sea = 0 and mass-derivative ∂ mq Z V (m q ). Combined they give estimates ofR V as listed in table 3. Its uncertainty is about one percent and below. For reasons described in the next section, we increase the uncertainty ofR V at β = 3.3 before probing various trial Padé functions in a weighted least square minimisation. The data is best described bȳ which we take as our preferred representation of b V + N fbV . With the accuracy achieved in the present paper, its uncertainty is not really of any practical relevance, except when we derive the uncertainty ofb V in section 2. (2.9) The functional form (2.8) is plotted in figure 3 together with the input data. We should remark that the improvement coefficient b V is known to one-loop order in perturbation theory [9,25], b PT V (g 2 0 ) = 1 + 0.0886(2)C F g 2 0 = 1 + 0.1181(3)g 2 0 . As the leading contribution ofb V is at O(g 4 0 ), we could constrain all Padé ansätze by b PT V (g 2 0 ). However, the accuracy of our data, especially at small couplings (β = 32.0, 16.0), reveals some tension to the one-loop estimate which is hard to account for in that case:p 1 −p 4 = 0.087(3). This could be due to the asymptotic nature of the perturbative series when compared to a nonperturbative result, or may expose a point we have not stressed explicitly yet. We have argued that the normalization conditions in use have hardly any systematic effects. On the other hand they explicitly depend on the non-perturbative clover-coefficient c sw [16] that renders the Wilson action O(a) improved. Although the Padé ansatz used for c sw (g 2 0 ) accounts for the correct one-loop coefficient, it undershoots it at couplings below g 2 0 ≈ 0.9 and approaches it only asymptotically. In that respect it is similar to what we observe but without data in the deep perturbative regime.
Of course, this would have gone unnoticed as any fit excluding our β = 32.0, 16.0 data can be easily constrained by the aforementioned one-loop b PT V . In the end, this discussion is of purely academic interest as no simulation of physical interest is ever performed in that region of couplings. But it is the reason why we do not insist in constraining our final result, eq. (2.8).

Valence quark mass dependence
Next we come to our determination of b V which is required in a partially quenched setup, e.g., with charmed valence quarks. For that purpose, we choose the ensemble closest to vanishing quark mass at each coupling, and for convenience measure f 1 and f V with hopping parameters κ val coinciding with the κ sea of other simulations listed in table 2. As a result, the measurements of Z V (m q ) for different valence quark masses are strongly correlated, leading to much more accurate estimates of R V . Even with higher accuracy, we do not observe mentionable differences from linearity in the range |am q | < 0.014 and as in the previous section stick to it in our final analysis. Before continuing we want to present a stringent test to explicitly validate the use of L/a = 8 lattices.  In figure 2 we show the Euclidean time dependence of R V at β = 4.0 where two volumes are at our disposal, T = L = 8a, 16a. The latter is available to us thanks to the work in ref. [18,26], and we refer to it for further details about that simulation. First, it confirms that, one or two lattice units away from the boundaries, there is no x 0dependence beyond statistical fluctuations, suggesting that O(a 2 ) effects in eq. (2.4) are small. Second, the agreement of R V (T /2) between both volumes shows that our estimate is volume independent. These observations hold within the statistical uncertainty which is at the per mille level or below. Monitoring the x 0 -dependence for all values of β, we remark that only at the coarsest lattice spacing (β = 3.3) we observe a +2σ deviation in both R V (T /2 + a) and R V (T /2 − a). To account for this still insignificant fluctuation, we conservatively double the uncertainty of R V (T /2).
We compile results of R V for two types of analysis in table 3. For reasons that become clear soon, we distinguish between a full (correlated) analysis and one neglecting correlations between measurements at different masses, which increases the uncertainty of R V by about an order of magnitude. Of course, both are in full agreement. In the course of simulations to produce ensembles at different sea quark masses we were guided by the statistical precision inR V , leading to a precision in the correlated analysis of R V that goes far beyond what is needed in future applications. While it allowed to probe relation (2.4) and the range of linearity more precisely than before, it is not straightforward to find a proper global approximation at that level of precision without increasing the number of data points.
For the sake of simplicity we thus prefer to employ data from the uncorrelated analysis on each ensemble, remarking once more that results at different couplings are statistically independent. At β = 4.0 we furthermore include the independent data point from L/a = 16. Not surprisingly, a (2.11) We add this interpolation and the entering data points to figure 3. One observes a significant difference toR V in the strong coupling region where (2+1)-flavor large volume simulations with the same bulk lattice action, produced by the Coordinated Lattice Simulations effort (CLS), exist [24,27,29]. Based on those simulations various estimates of mass-improvement coefficients have been published in ref. [30] for β ∈ {3.4, 3.46, 3.55, 3.7}. The uncertainty of the latter for b V is of the size of our estimated difference between sea and valence quark sector.

On sea quark mass-effects in the vector channel
With interpolating formulae for R V andR V we are finally able to check on the size of pure sea quark mass-effects that are cancelled at leading order by the countertermb V Tr [aM q ]. We defineb by eq. (2.8,2.10) and present it in figure 4.R V and R V are statistically independent, such that the uncertainty ofb V follows trivially from standard Gaussian error propagation. The value ofb V is different from zero in the non-perturbative region of large volume simulations [27], being about 0.07 at β = 3.4. For those simulations the unsubtracted trace of the  [23] this implies Tr [aM q ] ≈ 0.02. The resulting counterterm to cancel mass-dependent sea quark cutoff effects of a strange quark thus becomes an effect at the per mille level. For a dynamical heavy charm quark in the sea this effect would rise to the per cent level even without including higher order effects.
Given the similarity between b A , b P and b V , cf. [14], one could be lead to the assumption thatb A andb P are comparable in size to the result forb V presented above, but only a direct computation can confirm that.

Epilogue
In combination with a precise determination of the overall renormalization factor Z V and improvement coefficient c V along a line of constant physics [31][32][33], the determination of b V and b V presented in this paper completes the on-shell O(a) improvement of the nonsinglet local vector current with massive Wilson fermions that is consistent with massless renormalization schemes. Our results are independent from determinations of Z V and c V , and thanks to the applied improvement conditions we have reached very accurate estimates ofR V and R V that allowed for disentangling sea and valence quark effects with large significance for the first time. These estimates are particularly relevant to lattice phenomenology on CLS (2+1)-flavor ensembles with lattice spacings in the range 0.04 a/fm 0.09. Due to the smallness ofb V even at the non-perturbative level (a 0.1 fm) it seems likely that these effects are negligible in most applications of lattice QCD with dynamical light quarks. We have shown explicitly that higher order mass-dependent effects are practically absent, which is not true anymore for, e.g., dynamical charm quarks.
Acknowledgements. We thank M. Lüscher for a critical reading of the manuscript. P.F. acknowledges financial support by the Spanish MINECO's "Centro de Excelencia Severo Ochoa" Programme under grant SEV-2012-0249, as well as from the MCINN grant FPA2012-31686. The main part of this work has profited from the IFT computing infrastructure, especially the Hydra cluster, and some additional measurements were performed on the CERN cluster. We thank both IT departments for the provided resources and support.