Non-perturbative renormalization and O$(a)$-improvement of the non-singlet vector current with $N_{\mathrm{f}}=2+1$ Wilson fermions and tree-level Symanzik improved gauge action

In calculating hadronic contributions to precision observables for tests of the Standard Model in lattice QCD, the electromagnetic current plays a central role. Using a Wilson action with O($a$) improvement in QCD with $N_{\mathrm{f}}$ flavors, a counterterm must be added to the vector current in order for its on-shell matrix elements to be O($a$) improved. In addition, the local vector current, which has support on one lattice site, must be renormalized. At O($a$), the breaking of the SU($N_{\mathrm{f}}$) symmetry by the quark mass matrix leads to a mixing between the local currents of different quark flavors. We present a non-perturbative calculation of all the required improvement and renormalization constants needed for the local and the conserved electromagnetic current in QCD with $N_{\mathrm{f}}=2+1$ O($a$)-improved Wilson fermions and tree-level Symanzik improved gauge action, with the exception of one coefficient, which we show to be order $g_0^6$ in lattice perturbation theory. The method is based on the vector and axial Ward identities imposed at finite lattice spacing and in the chiral limit. We make use of lattice ensembles generated as part of the Coordinated Lattice Simulations (CLS) initiative.

In calculating hadronic contributions to precision observables for tests of the Standard Model in lattice QCD, the electromagnetic current plays a central role. Using a Wilson action with O(a) improvement in QCD with N f flavors, a counterterm must be added to the vector current in order for its on-shell matrix elements to be O(a) improved. In addition, the local vector current, which has support on one lattice site, must be renormalized. At O(a), the breaking of the SU(N f ) symmetry by the quark mass matrix leads to a mixing between the local currents of different quark flavors. We present a non-perturbative calculation of all the required improvement and renormalization constants needed for the local and the conserved electromagnetic current in QCD with N f = 2 + 1 O(a)-improved Wilson fermions and tree-level Symanzik improved gauge action, with the exception of one coefficient, which we show to be order g 6 0 in lattice perturbation theory. The method is based on the vector and axial Ward identities imposed at finite lattice spacing and in the chiral limit. We make use of lattice ensembles generated as part of the Coordinated Lattice Simulations (CLS) initiative.

I. INTRODUCTION
Precision tests of the Standard Model typically require reliable theory input from firstprinciples calculations. While the electroweak sector can be treated perturbatively, the virtual contributions of hadrons must be calculated from QCD non-perturbatively. Ab initio Monte-Carlo simulations of lattice QCD have provided a host of precision quantities for use in tests of the Standard Model [1]. Example of such hadronic quantities are the ratio of decay constants f K /f π , the MS quark masses, the running strong coupling constant α s (M Z ), and the anomalous magnetic moment of the muon, (g − 2) µ . For the latter, a major effort by several lattice collaborations worldwide is ongoing to calculate the hadronic vacuum polarization and the hadronic light-by-light contributions [2]. From the QCD point of view, these contributions amount to two-and four-point correlation functions of the electromagnetic current, to be integrated over with a weight function containing the characteristic scale of the muon mass.
In continuum QCD, the electromagnetic current is conserved and does not require renormalization. On the lattice, a finite renormalization can appear, depending on the details of the action and of the chosen discretization of the vector current. In particular, for Wilson fermions, the single O(a) on-shell improvement term to the action is known. Wilson fermions also have a 'point-split' vector current, whose support extends over two lattice sites in the direction of the current, which is exactly conserved at finite lattice spacing. This appealing property however does not guarantee that transverse correlation functions of the current have smaller discretization effects than those of the naive, "local" vector current with support on a single lattice site, which in the limit of massless quarks receives a finite renormalization factor Z V (g 2 0 ). Indeed, the improvement of the vector current -local or point-split -requires adding the divergence of the tensor current with a coefficient denoted c V , which counteracts the breaking of chiral symmetry by the Wilson action and suffices to remove all O(a) cutoff effects in on-shell correlation functions. This coefficient, whose value depends on the discretization of the current, has a finite value at tree-level of perturbation theory in the case of the point-split current, but vanishes for the local current.
On the other hand, for the local vector current, a mass dependence of the renormaliza-arXiv:1811.08209v1 [hep-lat] 20 Nov 2018 tion factor arises if O(a) discretization errors are to be removed. This mass-dependence is relevant in precision calculations, given the pattern of the physical up, down and strange quark masses. Concretely, given the electric charge matrix of the lightest three quark flavors, Q = diag(2/3, −1/3, −1/3), the electromagnetic current can be written as the linear combination V e.m.
where V a µ =ψγ µ λ a 2 ψ is the octet of vector currents, with λ a the Gell-Mann matrices. In isospinsymmetric QCD, the bare quark mass matrix can be decomposed as See Eq. (8) and below for our notation. The renormalization pattern of the local discretization of the two neutral octet combinations then reads [3], at O(a), with V 0 µ = 1 2ψ γ µ ψ the flavor-singlet current. Here V 3,I µ and V 8,I µ are understood to already contain the improvement term proportional to c V . All coefficients appearing in the two equations above are functions of the coupling g 0 . In this article, we present a non-perturbative determination of the renormalization factors Z V , b V and b V as well as of c V , while the coefficient f V will remain undetermined. As explained below, there are however good reasons to expect f V to be numerically very small [3]. The improvement coefficient c V is determined by imposing continuum chiral Ward identities, as proposed in quenched QCD in Ref. [4]. We follow the presentation of Ref. [3] for the full renormalization and improvement in large volumes with N f = 2 + 1 Wilson fermions. The mass-dependent renormalization with N f = 2 Wilson fermions has been computed in Ref. [5]. Note that the method of Ref. [6] only allows one to compute a linear combination of the improvement coefficients for the conserved and local currents, and is insufficient to provide a full improvement condition for either discretization.
Our main motivation for the present calculation is to determine the two-point function of the electromagnetic current with only O(a 2 ) discretization effects. This will in particular allow for a shorter continuum extrapolation of the leading hadronic contribution to the anomalous magnetic moment of the muon, and therefore a more cost-effective set of lattice QCD simulations. Given that phenomenologically, the π + π − channel, which is described by the timelike electromagnetic form factor of the pion, accounts for more than two thirds of the total hadronic contributions, it is very natural to impose the renormalization condition on the local vector current that the electric charge of the pion be unity at every lattice spacing. This is the main renormalization condition we will adopt to determine We begin by giving an overview of the required theory background, which allows us to define our notation. We present the setup for our numerical calculation in section III and the results in section IV. We finish with some concluding remarks in section V. Appendix A presents a determination of the improvement coefficient c A of the axial current, and appendix B contains some results on the employed correlation functions in lattice perturbation theory.

A. Definitions and notations
We use Euclidean Dirac matrices, {γ µ , γ ν } = 2δ µν . We consider initially the general case of N f flavours of quarks. Flavor indices will be denoted by latin letters i, j, . . . Let be the bare axial current and pseudoscalar density. The on-shell improved operators are given by where c A is an improvement coefficient. The average bare PCAC quark mass m ij of quark flavours i and j is defined through the relation We also defined the subtracted bare quark mass of flavour i, Often, the hopping parameter κ i ≡ (8 + 2am 0,i ) −1 is used to parametrize the bare quark mass m 0,i . The value κ cr ≡ (8 + 2am cr ) −1 of the hopping parameter is the value for which the PCAC mass, defined through Eq. (7), vanishes in the SU(N f )-symmetric theory. The bare quark mass matrix is defined as M 0 = diag(m 0,1 , . . . , m 0,N f ), and similarly for the subtracted bare quark mass matrix, M q = diag(m q,1 , . . . , m q,N f ). Finally, we also introduce the average quark masses Here we will be concerned with the improvement and renormalization of the vector current V (ij) µ on the lattice. Two discretizations are in common use, the local (l) and the point-split (c) lattice vector currents, Instead of the point-split vector current, we actually consider the symmetrized version (cs) which has the same properties under space-time reflections as the local vector current 1 [7]. This ensures that the same counterterms are present to remove O(a) artifacts with the local tensor current defined as and where we use the symmetric lattice derivative, Generically, the renormalization pattern of the quark bilinears, including O(a) mass-dependent effects, has been derived in Ref. [3]. For the vector current, and writing V µ as a flavor matrix, it reads 1 The authors would like to thank Stefan Sint for pointing out this fact is the modified bare coupling, which is in one-to-one correspondence with the lattice spacing, irrespective of the quark masses [8]. The symbol 'tr' refers to the trace over flavor indices and the λ is any element of the SU(N f ) Lie algebra. The improvement coefficients c V , b V , b V and f V are function of the bare coupling only; Z V has no anomalous dimension and does not depend on the renormalization scale. Given that the coefficient b g is so far only known perturbatively, it is worth noting the following. If one Taylor-expands the function Z V and only keeps terms up to O(a), the expression (15) is equivalent to replacing the argument of Z V by g 2 0 and then substituting b V by Therefore, the renormalization conditions we use for the vector current are only able to determine the combination b eff V . In a second step, using the perturbative estimate of b g , we obtain a value for b V . In the future, when a non-perturbative determination of b g becomes available, the value of b V can be updated.
In section II B, we describe the strategy used to determine the renormalization constant Z V and the improvement coefficients b V , b eff V and c V .

B. Vector Ward identities and determination of
We define an infinitesimal local vector transformation by where the matrix λ acts on flavour space. Using the path integral definition of an expectation value and noticing that the previous transformation is a change of integration variables with unit Jacobian, one obtains the following identity, where S is the Euclidean action and O is any operator. In fact, the equality holds on every single gauge-field configuration because only the fermionic part of the action is affected. For Wilson-Clover fermions, it leads to the well-known vector Ward identity [9] δO δα(y) where ∇ * µ φ(y) = (φ(y) − φ(y − aμ))/a is the backward lattice derivative in the µ-direction, λ denotes the matrix transpose of λ and V Working in components, we now consider the vector transformation δψ i (y) = +α(y)ψ i (y) , δψ i (y) = −α(y)ψ i (y) for one specific flavor i. Then, using O(x, z) = P (ji) (x)P (ij) (z) as a probe operator with i = j, one finds such that Eq. (20) reads Summing over the spatial vector y in Eq. (23), the spatial derivative does not contribute due to the use of periodic boundary conditions and only the time derivative remains. Therefore, the three-point correlation function a 3 y P (ji) (x)V c,(ii) 0 (y)P (ij) (z) , viewed as a function of y 0 , is a piecewise constant function with discrete steps of +1 at y 0 = z 0 and −1 at y 0 = x 0 . In particular, for x 0 > y 0 > z 0 , the ratio R defined by is unity such that the point-split vector current does not need any renormalization factor : On the other hand, the local vector current is not conserved on the lattice and needs to be renormalized.
In N f = 2 + 1 QCD with a quark mass matrix given by (2), by imposing either of the ratios to equal unity on the lattice at finite quark masses, one can determine the renormalization factor of the local isovector current , including the O(a) mass-dependent terms, as given explicitly in Eq. (3). We note that this renormalization condition does not require the knowledge of c V , and that the two choices for the "spectator quark" correspond to two different renormalization prescriptions. Using ensembles with different values of m q,1 = m q,2 and m q,3 = m q,s , each parameter can be determined independently. We remark that Z V , b V and b eff V could also be determined in the same way from the matrix element of since the flavor-singlet charge operator does not contribute on a kaon state. On the other hand, to obtain sensitivity to the coefficient f V , an external state with non-vanishing baryon number is required, for instance one may require the ratio to vanish. Without the vector current improvement term proportional to f V , R ∆ ++ would receive a contribution of order a from disconnected diagrams; the role of the coefficient f V , which multiplies the flavor-singlet vector current, under which the ∆ ++ baryon is charged, is to cancel this contribution. Therefore, the magnitude of f V is determined by the size of disconnected diagrams with the insertion of a single vector current. In perturbation theory, f V is therefore of order g 6 0 , because at least three gluons must be emitted from the quark loop 2 Figure 1: Illustration of the chiral Ward identity in the continuum and in the limit m 12 = 0. Continuous horizontal lines indicate that the operator is projected onto vanishing spatial momentum.

C. Axial Ward identities and determination of c V
Once the renormalization factor Z V and improvement coefficients b V and b eff V are known, the improvement coefficient c V can be determined by enforcing an axial Ward identity. In the continuum theory, the latter can be derived from the specific transformation As the operator to be chirally rotated, we choose O(y) = A (23) µ (y), and we have Choosing α(x) to be unity inside the slab t 1 < x 0 < t 2 and zero outside, the variation of the action (per unit α) is given by with t 1 < y 0 < t 2 . We perform the integral over the divergence of the axial current explicitly in the continuum using Gauss's theorem. With the additional constraint z 0 / ∈ [t 1 , t 2 ], we then obtain the following Ward identity valid in the continuum [4], and impose it to hold on the lattice, at fixed quark mass and bare lattice coupling, up to higher order corrections O(a 2 ). The Ward identity in the chiral limit is illustrated in Fig. 1. For each discretization of the vector current (α = l, cs), we define our estimator for the local vector current andẐ (13) V (g 2 0 , m av q , m q,13 ) = 1 for the conserved vector current. In Eq. (32) we use the symmetric lattice derivative, as in Eq. (14). Other choices would differ by an O(a) ambiguity in the definition of c V . The renormalized axial and pseudoscalar operators for i = j are respectively given by in terms of the improved operators defined in Eq. (6). The renormalized quark mass is defined through the relation such that the combination m R,12 P (12) R is insensitive to Z P , b P and b P . The renormalization factor Z A and the improvement parameters b A and b A have been determined non-perturbatively in Refs. [10][11][12].
In Eq. (32), the operator O ext can be either the vector k0 (z 0 , 0) and does not need to be O(a)-improved. In perturbation theory, the choice of the tensor operator is superior, since in the continuum, the right-hand side of the Ward identity vanishes in the chiral limit; on the lattice, the improvement term then involves the two-point function of the tensor current and its task is to cancel the vector-tensor correlation, which is O(a) and originates from chiral symmetry breaking at the cutoff scale. As we will see in the next section, in the non-perturbative QCD vacuum, both choices are equally well suited for separations between the operators of order of 0.5 fm, because the vector-tensor correlation is then non-vanishing even in the continuum.
There is one subtlety here. In Eq. (30), we sum over all timeslices in the range [t 1 , t 2 ] which implies the presence of a contact term for x 0 = y 0 . Therefore, on-shell O(a)-improvement is not sufficient to remove all O(a) contributions and the limit m 12 → 0 must be taken to remove this contact term. This is done by computing the effectiveĉ V for different light quark masses using Eq. (32) and then extrapolating to the chiral limit.
Finally, in appendix A we briefly describe a way to determine the improvement coefficient c A using an axial Ward identity. Our non-perturbative determination of c A , which we can compare to the literature [24], serves as a cross-check of our numerical setup.

D. Known perturbative results
The known perturbative results in QCD with N c colors and N f flavors of quarks are the following. The result b g = 0.012000(2) N f g 2 0 + O(g 4 0 ), independently of the pure gauge action, was obtained in [13]. For degenerate quarks, only the combination b V + N f b V appears, and the perturbative series for b V starts at order g 4 0 . For the tree-level improved Lüscher-Weisz action, the results are (C F = N 2 c −1 2Nc ) [14,15] The tree-level value of c cs V is 1 2 .

III. NUMERICAL SETUP
We use the N f = 2 + 1 CLS (Coordinated Lattice Simulations) lattice ensembles whose main parameters are given in Table. I. They have been produced using the openQCD code 3 of Ref. [17] using the Wilson-Clover action for fermions and the tree-level Symanzik-improved gauge action. The parameter c SW has been determined non-perturbatively in Ref. [18]. We consider four values of the bare coupling β = 3.40, 3.46, 3.55 and 3.70 which correspond to lattice spacings in the range 0.050 − 0.085 fm [16]. Ensembles using (anti)-periodic boundary conditions (PBC) and open-boundary conditions (OBC) in the time direction have been generated on three different chiral trajectories. Two trajectories with constant m av q and m q,s = m phys q,s can be used Table I: Parameters of the simulations: the bare coupling β = 6/g 2 0 , the lattice resolution, the hopping parameter κ, the pion mass m π and the PCAC mass.Ẑ   to extrapolate results to the physical limit with physical masses of the light and strange quarks. A third trajectory uses degenerate light and strange quarks with m q,l = m q,s . Concerning c V , it is enough to consider ensembles on a single chiral trajectory (e.g. m av q = constant). However, to determine the two improvement coefficients b V and b eff V , we have to consider at least two different chiral trajectories.
For the calculation of the renormalization factor Z V , we need to compute the following threepoint correlation function, projected on vanishing momentum and two-point correlation functions Correlation functions are calculated using U(1) stochastic sources with time dilution [19]. On each gauge configuration, we generate an ensemble of N s stochastic sources with support on a single timeslice and satisfying where each component is normalized to one, η a α (x) * For each stochastic source s with support in time-slice z 0 , we solve the Dirac equation and denote the solution vector Φ s i (x; z 0 ) = a 3 z S(x, z)η s i (z). Correlation functions are given by  with V = L 3 the spatial volume. We have used the γ 5 -hermiticity of the fermion propagator S(x, y) = γ 5 S(y, x) † γ 5 and Φ s ji is a sequential propagator given by Φ s ji (y; x 0 , z 0 ) = a 6 x,z γ 5 S j (y, x)γ 5 S i (x, z)η s (z). In practice, since the stochastic sources do not introduce a bias, the number of sources N s on each gauge configuration can be small. We choose N s = 12 such that the numerical cost would be the same if we used the usual point source method with a single source location.
To compute the correlation functions in Eqs. (30-31) we instead use point sources and the method of sequential propagators for the three-point correlation functions. A point source is first created on timeslice z 0 . Then, a sequential inversion is performed using the variation of the action between time slices t 1 and t 2 as a sequential source. We thereby have access to all y 0 values in the range [t 1 , t 2 ]. To increase statistics, we also average over equivalent polarizations k = 1, 2, 3. The values of t 1 , t 2 , z 0 and y 0 used in our simulations are summarized in Table II. We have computed the correlation functions entering Eq. (31) to leading order in lattice perturbation theory (see appendix B) in order to test our lattice QCD code.

IV. RESULTS
A. Results for Z V , b V and b eff V Away from the boundary, it is convenient to use the variables t = x 0 − z 0 and t 1 = y 0 − z 0 . For each ensemble, the value ofẐ V is estimated from the ratio of three to two-point correlation functions, defined through Eq. (24) with a local vector current. We choose j = l (spectator quark), which define our renormalization scheme. The ratio has the asymptotic behavior and is fitted to a constant in the plateau region where discretization effects are small. For ensembles with anti-periodic boundary conditions in time, we use C P V P (x 0 , y 0 ; z 0 ) → C P V P (x 0 , y 0 ; z 0 ) − C P V P (x 0 , y 0 + T ; z 0 ) to impose the vector Ward identity on each gauge configuration which can have a non-zero charge due to thermal fluctuations. Typical plots for the ensembles N200 and N300 are given in Fig. 2 and the results for all ensembles are summarized in Table I. In a second step, the renormalization constant Z V , and the improvement coefficients b V and b eff V at a given value of the bare coupling g 0 are obtained using the fit ansatẑ The ensembles included in the fit satisfy |am q,l | < 0.01 and am av q < 0.01 such that higher order corrections are expected to be small. The results for each value of β are given in Table III and the statistical error includes the error on κ cr . We note that the coefficient b V is significantly larger than the one-loop perturbative estimate given in Eq. (36b) and that b Finally, we perform a Padé fit to obtain the renormalization factor and the improvement coefficients as a function of the bare coupling g 2 0 . Our final results read which automatically reproduce the one-loop calculations and where the parameters and covariance matrices are given by Scaling for bV Scaling for bV Figure 4: Left panel: continuum limit of the ratio Z l V, sub /Z V where Z l V, sub was computed in [12] and where Z V refers to our own determination. Middle panel: b l V and b s V correspond to our determinations with the light or strange spectator quarks respectively and b F V is the value computed in [20]. Right panel: b F V is the value computed in [20]. We have defined the dimensionless parameter a = a/0.5 fm.
To ensure that O(a) ambiguities vanish smoothly toward the continuum limit, the renormalization of the vector current must be performed along a line of constant physics (LCP). Since the CLS ensembles have different volumes, we checked explicitly that the impact on the renormalization factor is extremely small. The observableẐ (12) V has been computed on two sets of ensembles (H105/N101 and H200/N202) generated using the same lattice parameters but with different volumes and the results quoted in Table I (6) respectively. The last results, where the source is close to the boundary, is slightly higher and might be affected by boundary effects. Those tests make us confident that with the procedure in Eq. (41) we indeed extract the matrix element in infinite volume.
As noted above, we could also choose j = s for the spectator quark and the values of the improvement coefficients b V and b eff V would differ by an O(a)-ambiguity. To study this effect, we have repeated the analysis with j = s and the results are given in Table III. We do not observe any difference for the renormalization constant Z V at our level of precision (both results should differ only by an O(a 2 ) ambiguity). For b V and b eff V , we observe variations by factors of at most 1.07 and 1.7 respectively. In Fig. 4, the continuum limit behavior of the ratio between the two results obtained using either a light or strange spectator quark is shown in blue. For b V , we observe the expected linear scaling with the lattice spacing. For b V , the ratio goes to one only if one includes higher order discretization effects, which appear to be sizeable.

Comparison of results with previous work
The renormalization factor Z V has been determined independently in [12] using the chirally rotated Schrödinger functional framework. In Fig. 4, we plot the ratio is extracted from [12] using the line of constant physics called L 1 and where the denominator corresponds to our own determination. This ratio goes rapidly to one in the continuum limit, even though the expected O(a 2 ) scaling is not observed. However, the maximum deviation, obtained at β = 3.40, is less than 1.6%. Empirically, the available data for the departure of the ratio from unity can be described by the sum of a linear and a quadratic term in the lattice spacing (not expected theoretically), or by the sum of a quadratic and a cubic term. The latter fit in fact describes the data slightly better, see Fig. 4. It also yields coefficients of reasonable size if one evaluates the lattice spacing say in units of 0.5 fm. The coefficients b V and b eff V have also been determined recently in Ref. [20] using a different setup, based on the QCD Schrödinger functional. A comparison with our results is given in Fig 6. For b V , we observe a deviation of about 5%, similar to the O(a) dependence on the spectator-quark estimated above. In Fig. 4, we show the continuum limit behavior of the ratio with our own results and we observe the expected linear scaling. However, for b eff V , the difference with the results quoted in Ref. [20] is significant, especially at large couplings g 2 0 . Again, as can be seen on Fig. 4, we do not observe a linear scaling in a for the ratio of the two determinations, and higher order corrections cannot be neglected. It suggests that this parameter suffers from a large ambiguity.
From a practical point of view, one should remember that a typical value of am av q is 0.005 on the β = 3.55 ensembles, so that with 3b eff V 0.16 even a 100% ambiguity on b V has an impact below the permille level. Conversely, it could be that O(a 2 ) effects compete with these terms, resulting in a substantial O(a) contamination in our determination of b eff V . For the physics applications discussed in the introduction, it is interesting to compare our values for the renormalization factorẐ (12) V of the isovector current to those one would obtain using the Z V values of [12] and the b V and b eff V values from [20]. We find that our direct estimates are always slightly lower, and that the relative difference depends almost only on g 2 0 : it is about 1.3% at β = 3.40, 0.8% at β = 3.46, 0.37% at β = 3.55 and 0.12% at β = 3.70. We conclude that these differences are of reasonable size, compatible with the expected a 2 (and higher order) ambiguity introduced by the choice of a specific renormalization condition.

B. Results for the improvement coefficient c V
For y 0 not too close from t 1 and t 2 , we can extract the value ofĉ V for each lattice ensemble. In practice, since we want to use a line of constant physics, we choose y 0 − z 0 = 0.77 fm and interpolate linearly between two neighbouring timeslices when necessary. Deviations from LCP due to the different sizes of the lattices used should be very mild since we are working in large volumes. The values of z 0 , t 1 and t 2 in Eqs. (30-31) as well as the values of y 0 used in the interpolation are quoted in Table II. Similar results are obtained using either the vector or the tensor operator as a probe operator in Eq. (31). In practice, we use the linear combination k0 (z 0 , 0) which helps to improve the statistical precision. For Z A , we use the results called Z l A,sub using the L 1 -LCP from [12]. For b A we used the published values in [10] and for b A we use values from [21] (in practice, we use b eff A which includes the b g -term for Z A , as in Eq. (17)). The results forĉ V for each ensemble are given in Table I and differ from c V by the presence of a contact term which vanishes in the chiral limit, as explained in Sec. II C. The improvement coefficient c V is obtained using a linear extrapolation in the light-quark PCAC mass m 12 at constant m av q and the results, for each value of the bare coupling, are summarized in Table III. As can be seen in Fig. 5, we observe a very mild chiral dependence. The error is dominated by the statistical precision of the correlation functions and the error on b A . The uncertainty on Z A and b A appears to have a negligible impact at our level of precision. Finally, we perform linear or quadratic fits in g 2 0 to determine c V as a function of the bare coupling. The    eff V we also compare our results with previous lattice determinations [10,20]. results for the local and the conserved vector currents read The parametrization is consistent with the perturbative predictions collected in Section II D. Our values for c l V are significantly smaller (in magnitude) than the preliminary values determined in [23] by applying a similar improvement condition in the Schrödinger functional. It could be due to a large O(a) ambiguity in the definition of c V . However, our values for Z V also differ by more that one standard deviation from the ones computed in [23]. Since Z V enters in the determination of c V , this could partly explain the disagreement. For example, if we use the preliminary results for Z V given in [23], our value would beĉ l V = −0.146 for H105 and c l V = −0.178 for N200. These observations highlight the need for a good control over Z V to determine precisely c l V . Since the point-split vector current is conserved, Z c V = 1, this issue is absent for the determination of c cs V .

V. CONCLUSION
We have determined non-perturbatively the renormalization constant and improvement coefficients of the local and point-split non-singlet vector currents with N f = 2 + 1 O(a)-improved Wilson quark action and the tree-level Symanzik-improved gauge action. Only one coefficient, f V , is missing but is also expected to be small, as it starts at O(g 6 0 ) in perturbation theory; in this regard, we note that for the two other mass-dependent improvement coefficients, b V and b V , the hierarchy expected from perturbation theory is indeed observed in our non-perturbative results. All these parameters are required for the full O(a)-improvement of the vector current  (24) and the reduction of discretization effects in lattice simulations. They are essential in the calculation of the hadronic vacuum polarization (HVP) contribution to the muon g − 2, where a precision below 1% is aimed at in the near future.
Full O(a)-improvement of the vector current requires to consider the renormalization factor Z V ( g 0 ) at the value of the renormalized coupling g 0 instead of the bare coupling g 0 . We have taken this difference into account by replacing the improvement coefficient b V by the effective parameter b eff V , thus avoiding the use of the unknown coefficient b g . We have obtained the renormalization factor and improvement coefficients by imposing vector and axial Ward identities at finite lattice spacing and bare quark masses on a set of large volume ensembles. Deviations from the line of constant physics in our renormalization scheme have been studied and shown to be small for Z V , b V and b V .
Our final results for the different β values used in CLS simulations are summarized in Table III. We also provide interpolating formulae through Eqs. (44) and (46). As a cross-check of our methods, we have recomputed the improvement coefficient c A and find good agreement with the results of Ref. [24], which employ an improvement condition set up in the Schrödinger functional.
Our calculation is the first non-perturbative determination of the improvement coefficients c α V with N f = 2 + 1 Wilson quarks for both the local and (the symmetrized version of) the point-split vector currents. The value for the local vector current is small and both values, for the local and point-split vector currents, are close to their perturbative values.
The comparison with the recent findings of Ref. [20] shows that a potentially-large O(a)ambiguity in b V remains, but that it should vanish smoothly in the approach to the continuum limit. For the vector current renormalization, we find important corrections to the expected asymptotic O(a 2 ) scaling for the difference between our results and the recent determination of Ref. [12]. However, we note the relative discrepancy is rather small, and smaller than that observed for two different normalization conditions for the axial current [11,12].
In the future, other improvement coefficients may be determined for the lattice action used here, thanks to the availability of an extensive set of CLS lattice ensembles. In particular, the N f = 2 + 1 hadronic contribution to the running of the weak mixing angle involves the flavorsinglet vector current, whose improvement coefficientc V is unknown. A method to determine the latter based on a chiral Ward identity was proposed in [3].