The renormalised $\mathrm{O}(a)$ improved vector current in three-flavour lattice QCD with Wilson quarks

We present the results of a non-perturbative determination of the improvement coefficient $c_\mathrm{V}$ and the renormalisation factor $Z_\mathrm{V}$, which define the renormalised vector current in three-flavour $\mathrm{O}(a)$ improved lattice QCD with Wilson quarks and tree-level Symanzik-improved gauge action. In case of the improvement coefficient, we consider both lattice descriptions of the vector current, the local as well as the conserved (i.e., point-split) one. Our improvement and normalisation conditions are based on massive chiral Ward identities and numerically evaluated in the Schr\"odinger functional setup, which allows to eliminate finite quark mass effects in a controlled way. In order to ensure a smooth dependence of the renormalisation constant and improvement coefficients on the bare gauge coupling, our computation proceeds along a line of constant physics, covering the typical range of lattice spacings $0.04\,\mathrm{fm}\lesssim a\lesssim 0.1\,\mathrm{fm}$ that is useful for phenomenological applications. Especially for the improvement coefficient of the local vector current, we report significant differences between the one-loop perturbative estimates and our non-perturbative results.


Introduction
When formulating lattice QCD with Wilson fermions, its violation of chiral symmetry is accompanied by the well-known characteristic that the discretisation errors of the theory receive contributions linear in the lattice spacing, a. Despite the conceptual advantages of Wilson fermions 1 , as a practical shortcoming this implies a typically slow convergence of physical quantities to the continuum limit. The obvious way to tame these linear discretisation errors -i.e., to approach the continuum limit by performing numerical simulations at a series of decreasing, very small lattice spacings -generally still represents a computationally hard task, particularly in the by now standard situation with three or four dynamical quark flavours.
According to the long-established recipe pioneered through Symanzik's continuum effective theory [1,2], lattice artifacts can be systematically removed order by order in a via adding appropriate higher-dimensional operators to the action and to the composite fields whose correlation functions are desired to extract observables of interest in lattice QCD. In fact, implementing this prescription by eliminating the leading, O(a) cutoff effects has great computational impact, because it accelerates the approach of physical quantities to the continuum in numerical simulations with Wilson's lattice QCD. If one restricts oneself to requiring improvement only for on-shell quantities [3], such as particle masses and matrix elements between physical states, the structure of the improved action for QCD and of the improved quark currents remains rather simple, see Ref. [4] for a detailed early account on the general theory. More precisely, while for a removal of O(a) cutoff effects in spectral quantities it suffices to improve the action by introducing the dimension-five Sheikholeslami-Wohlert (aka, "clover") term [5], with a non-perturbatively fixed coefficient c sw , the improvement of matrix elements of local operators necessitates the addition of higher-dimensional counterterms to them, together with the determination of associated improvement coefficients. The present work employs the lattice discretisation of QCD with N f = 3 mass-degenerate flavours of Wilson quarks and the tree-level Symanzik-improved gauge action [3]; the non-perturbatively tuned parameter c sw to improve this fermion action is available from Ref. [6].
In tandem with the issue of cutoff effects and the rate of convergence of the theory to the continuum limit, an unpleasant consequence of the explicit breaking of all chiral symmetries by the Wilson term in the action reflects in the Noether currents of chiral symmetry to be no longer protected against renormalisation. Therefore, the matrix elements of the axial and vector Noether currents between hadron states and the vacuum, which are parameterised by corresponding decay constants and constitute prime examples of physical quantities that can be accurately measured in Monte Carlo simulations, need to be correctly normalised. The main approach to fix the current normalisation constants consists in demanding continuum chiral symmetry relations to be fulfilled as normalisation conditions at finite lattice spacing [7,8]. Since then, it has proven itself in numerous settings to implement this strategy by exploiting chiral Ward identities, which follow from an infinitesimal chiral change of variables in the QCD path integral and, once evaluated numerically, are able to lead to a target precision that does not dominate the uncertainty of the (renormalised) currents' matrix elements.
In the following we focus on the non-perturbative improvement and renormalisation of the (iso-vector flavour non-singlet) vector current, since many interesting observables in lattice QCD are extracted from correlation functions involving this current in one of its most widely used discretisations, i.e., either the local vector current which has support on only one lattice site, or the conserved point-split one which extends over two neighbouring lattice points connected by a gauge link. Analogously to the renormalised O(a) improved axial vector current, which in the chiral limit of massless quarks is parameterised by an improvement coefficient c A and a normalisation factor Z A , the improvement and renormalisation pattern of these two lattice vector currents requires a single additive improvement term in each case, whose respective coefficient c V gives the size of its O(a) mixing with the (derivative of the) tensor current, multiplied by a finite (i.e., renormalisation scale independent) normalisation factor Z V . In the corresponding expressions, which for the vector channel are explicitly given in the next section, improvement coefficients and normalisations are functions of the bare coupling; c A and Z A , are non-perturbatively known for our action and phenomenologically relevant lattice spacings of about 0.1 fm and below from Refs. [9] and [10,11], respectively.
Our strategy to determine c V and Z V non-perturbatively orients on earlier investigations in the quenched approximation [8,12] and in two-flavour QCD [13] (Z V only). It adopts the Schrödinger functional formalism in conjunction with exploiting massive chiral Ward identities, which relate correlation functions of axial vector and vector currents. In the O(a) improved theory, and for vanishing quark mass, these identities can be written in a form that is valid up to error terms quadratic in a. Hence, the Schrödinger functional setup is particularly suited to numerically evaluate improvement and normalisation conditions inferred from Ward identities, because it enables to impose them in simulations close to or even in the chiral limit, in the very spirit of a mass independent scheme as it is adopted here. 2 Previous non-perturbative determinations of improvement coefficients and renormalisation factors, relying on the same Schrödinger functional gauge field configuration ensembles with three flavours of mass-degenerate sea quarks, include the aforementioned c A [9] and Z A [10], b-coefficients of additive quark mass terms as well as (ratios of) Z-factors for O(a) improved quark mass renormalisation [16,17], and the O(a) improvement of the tensor current [18]. With the present vector channel study to determine c V and Z V we complement these works. In the case of c V , for instance, our computation goes beyond the quenched one [12] also methodigically, as it more systematically includes, amongst a few more refinements on the technical level, the finite quark mass contributions to eventually eliminate them in a controlled way by an only slight extrapolation (or even interpolation) of the results to the zero mass limit in the unitary setup of equal sea and valence quark masses.
Other Ward identity calculations of Z V for almost the same coupling range of the three-flavour lattice QCD action at hand have recently been published in Refs. [11] and [19], but employ different frameworks. While the former makes use of the so-called 'chirally rotated' version of the Schrödinger functional [20,21], the latter directly operates on a subset of the large-volume N f = 2 + 1 gauge field configuration ensembles with open boundary conditions, generated within the joined effort by the CLS (Coordinated Lattice 2 Improvement conditions based on chiral Ward identities used in the past in pure gauge theory without Schrödinger functional boundary conditions can be found in Ref. [14]; see also [15] for its generalisation to lattice QCD with light dynamical quarks of non-degenerate masses. Simulations) cooperation of lattice QCD teams [22][23][24][25]. Moreover, [19] also addresses the O(a) improvement coefficients for the local and conserved vector currents discussed here (and defined below), as well as the related b-coefficients for the flavour non-singlet part of the local current that arise with the mass dependence of the renormalisation factor if O(a) discretisation errors are to be removed at physical values of the light quark masses.
In particular for the lattice prescription of the local current, estimates of the relative contribution of the improvement term evaluated with the perturbative one-loop value for the improvement coefficient c V known from [26] may suggest that the effect of improvement in correlation functions entering, e.g., decay constants would be small, and also the nonperturbative results found in Ref. [19] stay very close to this perturbative prediction, even in the region of stronger couplings. On the other hand, already the preliminary outcome of the early quenched calculation [12] in the Schrödinger functional setup indicated significant deviations of the non-perturbative estimates from perturbation theory towards the strongcoupling regime. And in fact, also in the present three-flavour investigation, we observe sizeable values for c V , which in the region g 2 0 1.55 are far from one-loop perturbation theory, whereas we recover the expected asymptotics in agreement with the perturbative behaviour when g 2 0 → 0. Another conceptual difference compared to Refs. [12,19] is that our improvement and renormalisation conditions based on chiral Ward identities are imposed along a line of constant physics (LCP). For this purpose, the spatial extent of the physical volume is fixed to L ≈ 1.2 fm, while T /L ≈ 3/2, such that all length scales in correlation functions are kept constant in physical units and only the lattice spacing a changes when g 0 is varied. 3 Since this is just the situation, to which Symanzik's local effective theory of cutoff effects can be applied, the resulting estimates of c V and Z V are expected to exhibit a smooth dependence on the bare coupling. Then it is obvious that any remaining intrinsic ambiguities, which could potentially emerge through another choice for the LCP or, equivalently, from a differently defined improvement or renormalisation condition, will thus asymptotically disappear towards the continuum limit at a rate of at least ∝ a in case of c-coefficients and ∝ a 2 in case of Z-factors if O(a) improvement is fully operative.
Since our range of bare gauge couplings covers lattice spacings 0.04 fm a 0.1 fm, which matches those of the aforementioned large-volume N f = 2 + 1 flavour QCD ensembles by CLS with the same lattice action, our results may be used for a variety of phenomenologically relevant applications involving the vector current, for which knowledge of c V and Z V is required. To name a few prominent examples, these not only comprise computations of vector meson decay constants, semi-leptonic decay form factors and the leading-order hadronic contributions to the muon's anomalous magnetic moment (a major part of which is accounted for by the timelike pion form factor) 4 , but also thermal correlators for the di-lepton production rate in the quark-gluon plasma. As another potential area of usage we mention the non-perturbative matching of Heavy Quark Effective Theory (HQET) to QCD in finite volume along the lines of Refs. [28,29], where the renormalised axial and vector currents enter on the QCD side of the matching calculation in order to extract the HQET parameters required for determining B-meson decay constants and semi-leptonic 3 For an earlier discussion and applications of this idea, see [27] and references therein. 4 In computing these hadronic contributions, the electromagnetic (quark) current plays a key rôle, which in QCD with the lightest three quark flavours is conveniently written as a linear combination of two neutral combinations from the octet of vector currents, V el.-mag.
form factors. The remainder of this paper is structured as follows. After a recapitulation of the theoretical framework in Section 2, in Section 3 we introduce the N f = 3 gauge field configurations underlying this work. Section 4 summarises the basic chiral Ward identities between quark currents, our analysis, and final results on the normalisation constant Z V of the vector current and compares them with those of Refs. [11,19]. As a byproduct, we also obtain the coefficient of the (valence) quark mass dependent piece of the renormalisation factor away from the chiral limit, b V , for one value of the lattice spacing to demonstrate consistency with the results of [19,30]. The Ward identity determination of the improvement coefficients for both, the local as well as the conserved (i.e., point-split) lattice discretisation of the continuum vector current is discussed in Section 5, finally yielding a continuous parameterisation of our non-perturbatively computed values as a function of the bare coupling. Since in the strong-coupling domain and especially for c V we find, by contrast to Ref. [19], sizable non-perturbative corrections to the one-loop formula, we present a few consistency checks that demonstrate higher-order ambiguities among these different determinations to vanish as g 2 0 → 0. Section 6 contains our conclusions; further data tables and definitions of Schrödinger functional correlation functions are deferred to appendices.
A preliminary account on the present work was reported in [31].

Definitions
To prepare the ground for the non-perturbative improvement and renormalisation of the vector current in Wilson's lattice QCD through numerical simulations, this section summarises the theoretical framework and a few relevant definitions.

Quark currents
We define the flavour non-singlet fermion bilinears associated with the vector current, tensor current, axial vector current, and the pseudoscalar density as where T a are the anti-hermitean generators of SU(N f ). In this work we consider three flavours of degenerate quarks with the bare subtracted quark mass 3) The corresponding current quark mass can be obtained via the PCAC relation as for any operator O a . In applications to hadronic physics one would also like to control the O(a) effects that according to Symanzik's programme applied to Wilson's lattice QCD are canceled by the counterterms to the currents. For the axial and vector currents, and close to the chiral limit, one essentially requires counterterms with coefficients c A and c V , respectively. More precisely, the improvement and renormalisation pattern of these currents reads where M q = diag(m q,1 , m q,2 , m q,3 ) = m q 1.
Alternatively one can define a conserved vector current, which is not local but split among neighboring lattice sites: In practice we employ a symmetrised version of the conserved vector current [32] defined as which under spacetime reflections transforms in the same way as the local vector current. 5 As the name suggests, the conserved vector current does not require renormalisation but local improvement via a counterterm parametrised by the coefficient cṼ For the renormalisation of the local vector current and the accompanying improvement coefficient there exist one-loop perturbative predictions for the tree-level Symanzik-improved gauge action [26,33]: For the improvement coefficient of the conserved vector current, cṼ, only its tree-level value 1 2 is known.

Schrödinger functional
For the numerical determination of the renormalisation and improvement factors of interest we evaluate the QCD path integral on a hypercubic space-time lattice with Schrödinger functional boundary conditions (i.e., periodic in space and Dirichlet in time). Since these introduce a gap into the spectrum of the Dirac operator, the use of the Schrödinger functional enables us to simulate with quark masses very close to the chiral point defined by am = 0. Another advantage can be gained by exploiting the property that the boundary time slices are invariant under global gauge transformations. The parity constraints at the boundary only allow for fermion bilinears that are parity odd. This essentially leaves us 5 In all situations where we consider the conserved vector current, we explicitly sum over all spatial lattice points. Due to periodic boundary conditions and zero momentum, the spatial components of the two variants of the conserved current are equivalent when the sum is taken: with two choices for the boundary sources, one corresponding to a pseudoscalar, the other to a vector operator: The explicit forms of the Schrödinger functional boundary-to-bulk and boundary-toboundary correlation functions relevant for constructing the estimators of improvement coefficients and renormalisation factors computed here are given in later sections and in Appendix B, respectively.

Gauge field ensembles and analysis aspects
The parameters of the gauge field configuration ensembles used in this study are summarised in Table 1. The ensembles are designed to lie on a line of constant physics (LCP) defined by a physical extent L ≈ 1.2 fm. As the original tuning was based on the perturbative (universal) two-loop beta-function as explained in [9], L is only constant up to O(g 2 0 ) effects. We explicitly confirmed by exploiting the results on the scale setting for our lattice action from [23] that these cutoff effects are consistent with O(a) (cf. discussion in Subsection 4.2). We thus expect our final results to be only affected at O(a 2 ) for renormalisation constants and O(a) for improvement coefficients, which is beyond the order we are interested in and treated as an ambiguity that extrapolates to zero in the continuum limit and is not accounted for in the final error budget. Further details on the numerical simulations to generate the gauge field ensembles in Table 1 can be found in refs. [9,10,[16][17][18]31], where they have already been the basis of previous work on non-perturbative improvement and renormalisation in three-flavour QCD by means of imposing Ward identities formulated in terms of Schrödinger functional correlation functions.
In numerical simulations by standard Monte Carlo algorithms (and if the physical volume is not so small that fluctuations of topology are a priori suppressed) it becomes increasingly difficult to sufficiently sample all topological sectors of QCD towards the continuum limit. This problem, often referred to as "topology freezing", implies that for the finer lattice spacings our simulations and the results extracted from them may suffer from critical slowing down of Q. To make sure that an incorrect sampling of the topological charge does not affect our data, we project expectation values to the topologically trivial sector by reweighting, as suggested in [35,36] and also applied in Refs. [9-11, 16-18, 31] in similar contexts before. Circumventing this algorithmic issue by the restriction to the trivial topological sector provides a viable and theoretically sound procedure, since it corresponds to imposing within the Q = 0 sector improvement and renormalization conditions, which derive from chiral flavour symmetries holding separately in each topological charge sector. For the definition to evaluate the topological charge, we follow the prescription adopted in Refs. [35,36]. While this projection to Q = 0 hence comes at the expense of typically somewhat larger statistical errors 6 and a slightly modified cutoff dependence, it is expected D1k4 24 3 × 35 3.81 0.137033 85008 0.696 9.6(1.9) Table 1: Summary of simulation parameters for the gauge configuration ensembles labeled by 'ID'. MDU denotes the total number of molecular dynamics units. P (Q 0 ) labels the percentage of configurations, for which the topological charge Q vanishes. τ exp indicates the exponential autocorrelation time used for the tail within the data analysis along [34]. It is estimated from the integrated autocorrelation time of the correlator F 1 on the longest Monte Carlo chain for each value of β. All configurations in a given ensembles are separated by 8 molecular dynamics units, except for ensembles A1k3 (4) and D1k4 (16). This is why τ exp for these ensembles differs from the others at the same value of β.
(and also confirmed in our data) to not induce a noticeable difference in the final results and their g 2 0 -dependence exhibiting a smooth approach to the continuum limit. For the error analysis we use (a python implementation of) the Γ-method [37], also accounting for effects of critical slowing down according to [34], incorporating the techniques of automatic differentiation as suggested in [38]. As estimate for the autocorrelation time τ exp of the slowest mode in a given Monte Carlo history (after Q = 0 projection), we employ the integrated autocorrelation time of the boundary-to-boundary correlation function F 1 defined in eq. (4.5), based on the assumption that, due to the parity invariance of the HMC algorithm, it is sufficient to look for the slowest mode in parity-even observables [34]. We determine τ exp for the longest Monte Carlo chain available for each value of β and conservatively consider this result as also valid for the remaining ensembles, see Table 1.
Only in a few cases, τ exp has to be adjusted by a factor stemming from different configuration spacings in molecular dynamics units (MDU).

Renormalisation of the vector current
We estimate the finite renormalisation of the vector current by demanding that the vector Ward identity in the O(a) improved theory is fulfilled up to O(a 2 ) corrections. In the present three-flavour case, this idea is numerically implemented within the finite-volume Schrödinger functional setup, similarly to previous calculations in the quenched approximation and in two-flavour QCD [8,13], while here in addition we also explore to fix the normalisation of the vector current via an alternative choice of boundary operators, yielding another non-perturbative estimate of Z V with different O(a 2 ) corrections. Other recently reported determinations of Z V in the three-flavour theory for the same discretisation and almost the same β-range employ the chirally rotated Schrödinger functional [11] and the large-volume gauge field ensembles with open boundary conditions generated by CLS [19].

Renormalisation condition and its evaluation
We consider infinitesimal vector rotations of the quark fields defined as For the present discussion we restrict ourselves to the SU(2) subgroup of SU (3). Based on the invariance of composite operators under these rotations, one can derive the vector Ward identity in an integrated form: where we limit the vector variation to all time slices where y 0 < x 0 and set z 0 > x 0 . We identify the two operators O and O with Schrödinger functional boundary interpolators, which only change their flavour structure but preserve their Dirac structure under small vector rotations with abc the totally anti-symmetric tensor of SU (2). Utilizing the particular types of Schrödinger functional boundary-to-boundary correlation functions defined as we arrive at two variants of the vector Ward identity, viz. from which Z V can readily be extracted. Since in the spirit of adopting a quark mass independent renormalisation scheme we ultimately work in the (unitary) limit of vanishing degenerate sea and valence quark masses, the terms proportional to am q and Tr[aM q ] could have already been dropped at this stage.
We evaluate these two different variants of the Ward identity for every time slice and identify a plateau as depicted in the left part of Figure 1. As for the plateau average, we decide to use the central third of the temporal extent of the lattice such that the plateau lengths are ensured to be kept approximately constant in physical units. Generally, the variant Z k V seems to be less sensitive to boundary affect and settles to a plateau at smaller distances from the temporal boundary. However, the statistical fluctuations are smaller by a factor ≈ 2 in the case of Z f V , which is why we prefer this definition for our final results. We explictly confirm below that settling on one of these two definitions of the renormalisation constant only amounts to ambiguities which vanish with a rate ∝ a 2 as expected in the O(a) improved theory.
As illustrated for representative examples in the right part of Figure 1, a linear function is perfectly justified to extract the renormalisation constants defined at zero quark mass 7 from the data at given bare coupling under investigation. Here, the use of Schrödinger functional boundary condition facilitates very small positive and even negative quark masses, which then sometimes promotes the chiral extrapolations to interpolations and thus renders the fits very reliable. For the fitting procedure we rely on the orthogonal distance regression algorithm [39], which takes into account errors in both the dependent and independent variables. This increases the statistical error of the value at the chiral point by up to 20% in comparison to the standard least-squares procedure. As argued above, we do not use any b-coefficients in the chiral extrapolation, because our simulations are sufficiently close to the chiral limit. The results for the individual ensembles and its chiral limits are summarised in Table 4.
In the left part of Figure 2, the g 2 0 -dependences of our two determinations are presented in comparison to the results of Refs. [11,19]. Our non-perturbative results can very well be parametrised by smooth interpolation formulae of a form asymptotically approaching the one-loop perturbative prediction [33] in the g 2 0 → 0 limit, where, as discussed before, we advocate in future applications to employ the estimator corresponding to x = f, with fitted coefficients and covariance matrix (4.9c) In the top right panel of Figure 2, we compare the scaling of the different determinations of the local vector current renormalisation constant Z V . The lattice spacings are extracted from a fit to the results of [23]. The deviation of the ratio of our two different determinations Z k V and Z f V from 1 can be very well described by a term quadratic in a with χ 2 /d.o.f. = 0.91, which confirms the expectation that the leading relative cutoff effects are dominantly of O(a 2 ). The relative cutoff effects between our determination Z f V and Z G V from [19] are best described by a fit to 1 + c · a with χ 2 /d.o.f. = 0.35. We can also model the data with a fit to 1 + c 1 · a 2 + c 2 · a 3 with χ 2 /d.o.f. = 1.14. A priori, a scaling behaviour with leading O(a) effects is unexpected, as both determinations are supposedly valid up to cutoff effects of O(a 2 ). For this reason we advocate the two-parameter fit, which is also depicted in Figure 2, and conclude that our determination agrees with the one of Ref. [19] up to the expected ambiguities.
Moreover, the relative scaling between the result from the chirally rotated Schrödinger functional Z χ V [11] and ours on Z f V shows O(a 3 ) effects of significant size in addition to the leading O(a 2 ) ones, which was also observed in [19] for their result. The deviation from Z χ V /Z f V = 1 is described best by both a quadratic and a cubic term with χ 2 /d.o.f. = 0.42, while other two-parameter fits are incompatible. This is perfectly in line with the theoretical expectation of relative cutoff effects of at least O(a 2 ), but it may have practical implication on the subsequent determination of the vector current's improvement coefficient which turns out to depend rather sensitively on the choice of the renormalisation constant. To corroborate this point, we take a closer look at the ratio (Z A /Z V ) 2 (i.e., the relevant ratio entering the expression for c V , as we will see in the next section) for different variants   of Z V in the bottom right panel of Figure 2. While the ratio of the two determinations based on the same method seems to deviate from unity by effects dominated by O(a), its "mixed" ratios appear to exhibit additional ambiguities of O(a 2 ) that dominate for the coarser lattice spacings of interest.

Check of the LCP condition
Finally, we want to assess the systematic error that could be induced in our results by a violation of the LCP condition. Before discussing its quantitative outcome, let us recall the idea behind it: Imposing improvement and renormalisation conditions based on chiral Ward identities implies chiral symmetry to be fully recovered only in the continuum limit. Consequently, the choice of these conditions matters at the cutoff level, i.e., at a fixed value of the lattice spacing, the numerical results may differ -sometimes even considerablybetween any two such choices. However, when adhering to the prescription of first defining an LCP (by fixing all dimensionful parameters in terms of a physical scale) along which the continuum limit is taken, these intrinsic differences should not be interpreted as a systematic error. Rather, since this LCP approach yields functions c(g 2 0 ) and Z(g 2 0 ) as the lattice spacing is varied and, accordingly, different functions for another choice of the LCP, the difference between any two of these choices will be, within errors, a smooth function of g 2 0 that vanishes asymptotically at least ∝ a or ∝ a 2 if O(a) improvement is implemented. In other words, following an LCP ensures that cutoff effects obey a smooth functional dependence on g 2 0 and the particular choice of LCP becomes irrelevant in the continuum limit. Hence, from this perspective, the relevant systematic error is determined by the precision, to which a chosen LCP is realised and can be followed.
As the LCP used in this work was originally defined based on the perturbative two-loop β-function, we expect the deviation from the intended value of L in physical units to be proportional to the lattice spacing a. Information, on how accurately the chosen LCP condition is met for our simulation parameters, can be gathered from the left panel of Figure 3. The deviation can indeed be described by a function linear in the lattice spacing a. Based on this result we expect any effects of this deviation on our quantities of interest  to be of the next-to-leading order in a (i.e., O(a 2 ) for the renormalisation of the vector current in the O(a) improved theory).
To get a feeling for the size of these effects we explicitly investigate the volume dependence of Z V at a fixed value of β = 3.676. For this purpose we have generated additional gauge ensembles, summarised in Table 2, which coincide with the ensembles labeled C in Table 1 except for the lattice geometry. We evaluate Z f V at the chiral point on these ensembles as described in the previous subsection. Our estimates on the impact of and possible systematic uncertainties due to deviations from the chosen LCP are displayed in the right panel of Figure 3. Apart from the the result at L/a = 8, the data can be well described by a constant plus a term proportional to (a/L) 4 . The vertical dashed lines indicate the span of volumes which are part of the LCP. The empirical finding is that volumes of intermediate extents ∼ 1.2 fm, employed in this study, correspond to a regime where effects owing to a not exactly fixed physical volume 8 are mild but still of a similar magnitude as the statistical uncertainties achieved here. In the spirit of the previous line of reasoning, however, we decide to not include this uncertainty into our final results, because the deviations from the LCP are of O(a 2 ) or higher and therefore beyond the level that would have to be considered as a contribution to the systematic error.

Mass dependent O(a) improvement
For the sake of completeness we would like to briefly address the quark mass dependent pieces in the renormalisation patterns of the vector current according to eqs. (4.7) and  (4.7). To achieve full O(a) improvement of the local current at finite quark mass, one not only requires (besides c V ) the renormalisation constant Z V in the chiral limit, but also the coefficientsb V and b V which parametrise the sea and valence quark mass dependence, respectively. Both coefficients are already non-perturbatively known for our choice of action from Refs. [19,30]. Moreover, as the mass range width of our gauge ensembles is not ideal for a detailed study of these mass dependent effects for all g 2 0 , we restrict ourselves to an accurate determination of b V at only one value of the bare coupling, corresponding to β = 3.414, which we compare to the previous results. The coefficient b V is readily accessible through evaluating the logarithm of the renormalisation condition for Z V at different valence and vanishing sea quark masses. The slope with respect to the bare quark mass is then proportional to the b-coefficient of interest. Thereby, from the above renormalisation conditions (4.7) and (4.8), one arrives at two independent estimates of b V in terms of ratios of correlation functions by mapping out the valence quark mass dependence at vanishing sea quark mass: For the numerics we choose ensemble E1k2 (β = 3.414), which fulfills the condition of vanishing sea quark masses within statistical precision, and calculate the relevant ratios of correlation functions at five different, equally spaced valence quark masses, as shown in the left part of Figure 4. In the right part of Figure 4 we compare our two determinations with the results from the literature. Our b f V is very closely consistent with the results of [30], which emerge from an almost identical computational setup but employ a smaller physical volume. This supports the claim made in [30] that the volume dependence of the b-coefficients is not statistically significant. Although our result from the other estimator, b k V , has a considerably larger statistical uncertainty, it is perfectly compatible with the result of [19], in particular when keeping in mind that intrinsic O(a) ambiguities among the b-coefficients from different improvement (and realisations of) constant physics conditions must be expected in any case.

O(a) improvement of the vector current
The method for the determination of the vector current improvement coefficient based on chiral Ward identities has first appeared and been tested in [12] in the quenched case.
Recently, it has also been applied to QCD with N f = 2+1 dynamical fermions [19], using the same O(a) improved lattice action as we do, but implementing it directly on large-volume ensembles generated by the CLS effort. In this section we will briefly recapitulate the derivation of the improvement condition within our intermediate-volume approach realised by the Schrödinger functional setup and then propose a slight variation that helps to tame the large higher-order ambiguities in (Z A /Z V ) 2 , anticipated in the previous section, which now become influential.

Improvement condition
We now consider infinitesimal axial vector variations of the quark fields The improvement condition studied in the following is based on the invariance of the expectation value of composite operators under these rotations which is basis for the axial vector Ward identity presented here in its integrated form where we restrict the axial variation to the space-time region between t 1 and t 2 and impose z 0 < t 1 < y 0 < t 2 . When we identify O c with Q c k and O b with A b k , which changes flavour as well as Dirac structure under small axial vector rotations we arrive at .
The finite renormalisation constants Z A and Z V and the O(a) improvement terms proportional to c A and c V , respectively, assure that the Ward identity is valid up to quadratic corrections in the lattice spacing. (Explicit expressions of the Schrödinger functional boundary-to-bulk correlation functions involved are summarised in Appendix B.) This Ward identity can be rearranged to yield a first improvement condition for c V , viz., .
For the moment we state the functional (resp. parametric) dependence on the Euclidean times x 0 , t 1 and t 2 as well as on the renormalisation and improvement factors Z A , Z V and c A explicitly. The superscript 'I' labels improved correlation functions as detailed in Appendix B. In this improvement condition, c V arises as the difference of two potentially large terms, where only one of the two is sensitive to a change in Z A . Therefore, even a sub-permille uncertainty on Z A can easily propagate into a relative error well above five percent on the final error and may hence become the dominant contribution to the overall uncertainty on c V at the end. In order to suppress the relative higher-order effects observed in the determination [11] of the normalisation constants Z A and Z V via the chirally rotated Schrödinger functional, we propose an alternative improvement condition which exploits the previously discussed vector Ward identity to expand the fraction in the first term on the r.h.s. of eq. (5.6). This gives another estimator for the improvement coefficient c V , labeled by c V,alt , However, eq. (5.8) benefits from the fact that the ratio Z 2 A /Z 2 V , being directly deducible from Ref. [11], has balanced powers of the Z-factors (in contrast to Z 2 A /Z V entering (5.6)), from which one may expect an approximate cancellation of the dominating higher-order ambiguities as suggested by the dashed line in the bottom part of Figure 2 and thus also a reduced impact of these systematics on our result.
For the improvement coefficient of the point-split vector current, the two variants of the improvement condition look slightly different. For the original improvement condition we simply substitute k V by kṼ and use the fact that ZṼ = 1:   To arrive at an alternative improvement condition, where in a similar spirit as before the ratio Z 2 A /Z 2 V enters, we need to multiply the vector Ward identity (5.7) to the first term on the r.h.s. twice. This yields: (5.10)

Analysis details
We evaluate the improvement conditions for the local and the conserved (i.e., point-split) vector current with two different operator positions t 1 = T − t 2 ∈ {T /4, T /3} and average over the two central values of x 0 . For the axial vector improvement coefficient c A we employ the result of [9], while for the vector current normalisation constant we try both, the results Z f V from our determination reported here and the values Z χ V obtained within the chirally rotated Schrödinger functional setup [11]. For the axial vector current renormalisation constant we only rely on the result of [11] labeled Z χ A , because the accuracy reached for this quantity in the standard Schrödinger functional setup [10] is not sufficient for a satisfactory precision of the vector current improvement constants and would cap the final absolute errors above ≈ 0.1. Hence, from hereon we simplify the notation specifying our particular choices to define the two improvement coefficients for each of the two vector currents to Note that in the case of the standard improvement condition for the conserved (i.e., point-split) current, eq. (5.9), no Z V -dependence arises.
x 0 -dependences of the improvement conditions on an exemplary ensemble are presented in the left part of Figure 5. The fact that the local estimator does not develop a clear plateau but rather slightly descends towards larger values of x 0 was also seen in the similar study [19] as well as in the computation of the tensor current improvement constant [18]. We attribute this to the smaller Euclidean time distances of the operators contributing to the Ward identities in comparison to the PCAC relation that usually exhibits more distinct plateaux. Note, however, that by virtue of basing the extraction of improvement coefficients (and normalisation constants) on Ward identities formulated on the operator level and combining it with the constant physics idea, there is a natural freedom in choosing x 0 and t 1 when actually evaluating the given improvement condition in a numerical calculation, as long as all physical scales are kept fixed among the different ensembles. Therefore, even if the x 0 -dependence of the estimators does not exhibit a plateau in the central region, taking an average over central values of x 0 as implied by eq. (5.11) is perfectly legitimate. In case of the operator position t 1 = T /3, the plateau-like behaviour spans a shorter time separation as expected. For the reason of larger separation from the boundary and because of the slightly smaller statistical uncertainties, we prefer t 1 = T /3 over t 1 = T /4 and thus regard this choice as part of the definition of our estimators for c V ; the latter choice (i.e., t 1 = T /4) is only used for scaling analyses. In the right part of Figure 5 the chiral extrapolation for one specific improvement condition, c , is displayed. We compare the extrapolation with and without the explicit mass term (see, e.g., eq. (5.5)) to illustrate that accounting for this term in the evaluation of the Ward identity leads to an almost flat quark mass dependence and a thereby very stable extrapolation. The remaining slope could in principle be compensated by a factor involving mass dependent pieces proportional to b A ,b A , b V andb V . As a sufficiently precise non-perturbative determination of the mass dependent improvement coefficients for the axial vector current is not available for our setup, we decide to omit this factor which has no effect on the result in the chiral limit. The results for the individual ensembles and chiral extrapolations are summarised in Tables 5  and 6 in Appendix A.
In the left part of Figure 6 we compare three different determinations of c V from this work with the results of Gérardin et al. [19] and the one-loop prediction of [26]. While the results obtained by Gérardin et al. are compatible with perturbation theory, all three of our determinations strongly deviate from this and only come close to the expected one-loop prediction for c V as g 0 → 0. In particular, these three different variants agree within errors for the two smallest lattice spacings and show larger (albeit monotonic) spreads for the coarser ones. The two determinations c , which are sensitive to the higher-order ambiguities in Z χ A , show much milder scaling with g 2 0 compared to c with balanced powers of Z-factors, where the leading higher-order effects are supposed to be almost canceled. This suggests that these higher-order effects counteract the leading O(a) effects and imply the observed flat, somewhat non-uniform scaling behaviour with    [19] and one-loop perturbation theory [26]. Right: g 2 0 -dependence of different determinations of cṼ in comparison to the results of [19] and the tree-level value. g 2 0 , both present in the result of [19] and our determination c . For these reasons we finally decide to advocate c as our preferred determination and expect it to be least affected by ambiguities in the ratio Z A /Z V .
A similar observation can also be made for the g 2 0 -dependence of the improvement coefficient, cṼ, of the conserved vector current displayed in the right part of Figure 6. Since the standard improvement condition does not depend on Z V , we display only two variants here. As in the case of the local current, our results deviate much more from the tree-level formula in comparison to the result of [19]. Again we observe a rather mild and less smooth dependence on the (squared) bare coupling for the result of Gérardin et al. and our determination, c {T /3} V , which we attribute to the influence of higher-order effects in Z χ A . Therefore, also in this case we decide to prefer our alternative determination c over the former, for which we can well assume the dominant higher-order ambiguities to largely cancel.

Final results
The left part of Figure 7 collects our preferred determinations of the improvement coefficients of the local and conserved vector currents, respectively, in comparison to the perturbative predictions and the results of [19]. For brevity we refer to these as c V and cṼ. The figure also displays the corresponding continuous interpolations of these final non-perturbative results in terms of g 2 0 . The formulae are motivated by the leading term in the perturbative relation between the lattice spacing and the β-function, with universal coefficient b 0 = 9/(4π) 2 for N f = 3, constrained by the one-loop [26] (tree-level) perturbative prediction for g 2 0 → 0 in case of c V (cṼ). The g 2 0 -dependence is then best described by a parametrisation of the    [19] labeled by the superscript 'G'. The shaded points at g 2 0 = 0.75, i.e., lying deeply in the perturbative domain, are not included in the fits, as they do not lie on the LCP, but provide evidence that the improvement conditions approach the perturbative results towards the continuum limit as expected. Top right: Continuum limit of the difference of various determinations of cṼ with our preferred variant. Bottom right: Continuum limit of the difference of various determinations of c V with our preferred variant. In both figures on the right, the dashed lines describe the functional forms that describe the lattice spacing dependence best, as detailed in the text.
following form: (5.12d) The couplings used in our calculations largely overlap with the bare couplings of the N f = 2 + 1 flavour lattice QCD ensembles generated by the CLS consortium [22][23][24][25]. These large-volume gauge field configuration ensembles with open and periodic boundary conditions emerge from large-scale numerical simulations with the same lattice action and at almost the same range of lattice spacings as in the present work. Since they are designed for a variety of phenomenological lattice QCD applications, for the sake of completeness and for prospective use in future calculations involving the renormalised and O(a) improved vector currents, we also interpolate (respectively, in the case of β = 3.85, slightly extrapolate) our results to the corresponding values of β = 6/g 2 0 employed there. The corresponding estimates for Z V (resorting to the parameterisation (4.9) from Section 4), c V and cṼ are given in Table 3.  Table 3: Z V , c V and cṼ results for the inverse gauge couplings β = 6/g 2 0 of the N f = 2 + 1 CLS simulations [22][23][24][25]. The errors are the statistical uncertainties propagated from the interpolation formulae (4.9) and (5.12), except for β = 3.85, which lies slightly outside of the coupling range covered in the present investigation, where we add 50% of the size of the statistical error in quadrature as a systematic uncertainty accounting for this.

Consistency checks and scaling tests
We close this section with a discussion of a few consistency checks and scaling tests to demonstrate the vanishing of possible O(a) ambiguities.
First of all, as for the same check of volume dependence as for Z V in Subsection 4.2, we observe for the improvement coefficients a qualitatively very similar picture that (on the level of larger statistical errors) even points to a statistically insignificant L/a-dependence. Hence, the general conclusion drawn there also applies here. Furthermore, as an additional validation of our parameterising fit ansaetze in eqs. (5.12), we evaluate the improvement conditions in the weak-coupling regime, i.e., on a lattice of size 16 3 × 24 with g 2 0 = 0.75, and at the chiral point. Z A and Z f V are extracted along the lines of [10] and Section 4 of this this work, respectively, while c A is set to the one-loop prediction of [40]. As already mentioned before, the resulting numbers are also displayed in the left part of Figure 7. It is important to note, though, that this lattice has a distinctly smaller volume and does not lie on the LCP underlying the computations in the present work. Hence, although we do not expect these data points to perfectly comply with our interpolation formulae, we nevertheless see them as qualitatively very convincing support that our results in the regime of non-perturbative couplings, obtained via the Ward identity method, approach the perturbative predictions towards smaller bare couplings and that an ansatz monotonically connecting these two regimes is perfectly justified.
Next, in the top right part of Figure 7, we depict the outcome of a scaling analysis of different determinations of the conserved vector current improvement coefficient, cṼ. In the following discussion we exclude the data points at the coarsest lattice spacing, corresponding to β = 3.3 and outside of the CLS region, in order to describe the data with a reasonable functional form. The deviation of the difference of the two determinations c from zero can be described by a sum of terms quadratic and cubic in the lattice spacing. This is expected, as both are based on the same correlation functions and only differ by how the first term in the improvement condition is renormalised. The difference of the two determinations c exhibits much smaller corrections to zero that can reliably be fitted by a purely quadratic term in the lattice spacing. Since the operator placement t 1 of A 0 in the bulk is the only difference between the two, this is in line with our expectation that no dominant higher-order ambiguity is present here. The scaling of the difference of the result by Gérardin et al. [19] and our preferred determination is also well described by a quadratic ansatz in the lattice spacing of the form 0 + c · a 2 . In essence, intrinsic ambiguities among different definitions of cṼ-estimators, and even among completely independent determinations, vanish with a rate of at least O(a) towards the continuum limit as theoretically expected.
Similarly, the scaling analysis of different calculations of the local vector current improvement coefficient, c V , is illustrated in the bottom right part of Figure 7. Also here we exclude the data points at the coarsest lattice spacing. The deviation of the difference of our two determinations c from zero can again be fitted well by a sum of terms quadratic and cubic in the lattice spacing. The difference of c which only differ by the bulk position of the A 0 operators in the k I AA -andk PA -correlators, can be described by a term quadratic in the lattice spacing. For the difference of the result from [19] and our preferred determination, we find leading corrections that are linear in the lattice spacing, also in line with the theoretical expectation.

Conclusions
In the present investigation we have determined the O(a) improvement coefficients of the local and conserved (i.e., point-split) versions of the lattice QCD vector current, as well as the Z-factor fixing its normalisation, for the theory discretised à la Wilson. Our non-perturbative calculations are performed in the framework of lattice QCD with N f = 3 flavours of mass-degenerate, non-perturbatively O(a) improved [6] Wilson-clover sea quarks and tree-level Symanzik-improved [3] gluons. The underlying methodology is derived from massive chiral Ward identities, which are formulated and evaluated through numerical simulations in small, fixed physical volume of spatial extent L ≈ 1.2 fm and T /L ≈ 3/2, with Schrödinger functional boundary conditions and an unitary setup of valence quark masses equal to those of the sea quarks. In this way we are able to approach the chiral and continuum limits, while staying on a line of constant physics (LCP) when the Ward identities are imposed to be restored up to O(a 2 ) at finite lattice spacing. Since the Ward identities are valid as operator relations in all sectors with fixed topology, the correlation functions involved have been projected to the sector of zero topological charge in order to ensure a reliable statistical error estimation from Monte Carlo simulation data facing topology freezing.
Our main findings are the parameters of the interpolation formulae (4.9) and (5.12) for Z V and c V , cṼ, respectively, which hold for the range of bare couplings 1.55 g 2 0 1.85, i.e., lattice spacings 0.042 fm a 0.105 fm. It also covers the g 2 0 -values of the large-volume N f = 2 + 1 simulations by CLS [22][23][24][25], and the corresponding results at these couplings are collected in Table 3. 9 Note that in this strong-coupling region and for the improvement coefficients we observe sizable non-perturbative corrections to the one-loop (resp. tree-level) prediction, by contrast to the outcome of the study [19] employing a different improvement condition. Since the uncertainties on c V and cṼ receive a significant contribution from the (independent) statistical error on the ratio Z 2 A /Z V and may suffer from systematics due to unbalanced higher-order ambiguities between numerator and denominator, we have advocated a prescription where this ratio is replaced by Z 2 A /Z 2 V , very accurately know from Ref. [11], such that these effects are largely reduced. At this point, however, we would like to highlight again that, in the very spirit of the generic LCP idea, the exact choice of c-coefficients and Z-factors from our results is not decisive from the theoretical perspective as long as one stays within the set of numbers for a given definition and does not switch between different definitions when changing the bare coupling. Otherwise cutoff effects are no longer guaranteed to vanish smoothly at a rate ∝ a 2 .
Thanks to the constant physics condition, our results for the improvement coefficients and the normalisation factor exhibit a smooth dependence on the (squared) bare gauge coupling, monotonically linking the non-perturbative domain, where they typically are relevant for further applications in various contexts of hadronic physics, with the perturbative regime. Moreover, to demonstrate the negligibility of possibly uncontrolled systematic effects, we have demonstrated explicitly that variations in the actual numerical implementation of the Ward identity based improvement and normalisation conditions only give rise to ambiguities of expected higher order (viz., at least O(a) for c V and cṼ, and O(a 2 ) for Z V ) that uniformly disappear as g 2 0 → 0. We also find overall consistency of our results with those obtained in Refs. [11,19] by similar methods, though grounded on different ensembles and (in case of [19]) not strictly obeying the condition of constant physics. While for Z V the respective g 2 0 -dependences fall quite close to each other, c V and cṼ clearly disagree quite substantially from [19] at couplings over the common range of finite lattice spacings 0.04 fm a 0.1 fm, even though these differences are perfectly in line with O a (n≥1) ambiguities that as a consequence of completely independent improvement and renormalisation conditions are inherently unavoidable but vanish towards weaker couplings. In particular, together with the result for c V in Ref. [19], the outcome of our determination indicates that its one-loop perturbative prediction is compatible with both results, but up to significantly different higher-order cutoff effects. Therefore, it remains to be seen which set of estimates will in practice lead to a more pronounced impact of improvement in the scaling behaviour of quantities involving the vector current; this could, for instance, be quantitatively explored in a dedicated precision intermediate-volume scaling test similar to [41] or in future phenomenological studies. 10 As for the latter, interesting applications of the renormalised O(a) improved vector current and its matrix elements include computations of vector meson decay constants (such as, in the charm sector, the D * or J/ψ decay constants), semi-leptonic decay form factors, the leading-order hadronic vacuum polarisation and light-by-light contributions to the muon's anomalous magnetic moment, as well as thermal correlators related to the di-lepton production rate in the quark-gluon plasma.  Summary of results for the PCAC quark mass am and the two determinations of the vector renormalization constant Z f V and Z k V on the individual ensembles and in the chiral limit. The errors for the individual ensembles are statistical, the ones in the chiral limit follow from the orthogonal distance regression procedure of Ref. [39]. Our preferred determination Z f V is emphasized in italic. Table 5: Summary of results for am and different determinations of the local vector current improvement coefficient c V on the individual ensembles and in the chiral limit. The errors for the individual ensembles are statistical, the ones in the chiral limit follow from the orthogonal distance regression procedure of Ref. [39]. Our preferred determination c  Summary of results for am and different determinations of the conserved vector current improvement coefficient cṼ on the individual ensembles and in the chiral limit. The errors for the individual ensembles are statistical, the ones in the chiral limit follow from the orthogonal distance regression procedure of Ref. [39]. Our preferred determination c

ID
is emphasized in italic.

B Schrödinger functional correlation functions
The Appendix summarises the definitions of the Schrödinger functional correlation functions employed in this work.