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

###### Abstract

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_{\rm 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_{\rm 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_{\rm 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.

^{†}

^{†}preprint: MITP/18-106

## I Introduction

Precision tests of the Standard Model typically require reliable theory input from first-principles 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 Aoki:2016frl (). Example of such hadronic quantities are the ratio of decay constants f_{K}/f_{\pi}, the \overline{\rm MS} quark masses, the running strong coupling constant \alpha_{s}(M_{Z}), and the anomalous magnetic moment of the muon, (g-2)_{\mu}. 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 Meyer:2018til (). 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_{\rm V}(g_{0}^{2}). Indeed, the improvement of the vector current – local or point-split – requires adding the divergence of the tensor current with a coefficient denoted c_{\rm 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 renormalization 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, {\cal Q}={\rm diag}(2/3,-1/3,-1/3), the electromagnetic current can be written as the linear combination

V_{\mu}^{\rm e.m.}=V_{\mu}^{3}+\frac{1}{\sqrt{3}}V_{\mu}^{8}\,, | (1) |

where V_{\mu}^{a}=\bar{\psi}\gamma_{\mu}\frac{\lambda^{a}}{2}\psi is the octet of vector currents, with \lambda^{a} the Gell-Mann matrices. In isospin-symmetric QCD, the bare quark mass matrix can be decomposed as

M_{\rm q}=m_{\rm q}^{\rm av}+\frac{1}{\sqrt{3}}(m_{{\rm q},l}-m_{{\rm q},s})% \lambda^{8}\,. | (2) |

See Eq. (8) and below for our notation. The renormalization pattern of the local discretization of the two neutral octet combinations then reads Bhattacharya:2005rb (), at O(a),

\displaystyle V_{\mu,R}^{3} | \displaystyle= | \displaystyle Z_{\rm V}\;(1+3\overline{b}_{\rm V}\;am_{\rm q}^{\rm av}+b_{\rm V% }\;am_{{\rm q},l})\;V_{\mu}^{3,I}\,, | (3) | ||

\displaystyle V_{\mu,R}^{8} | \displaystyle= | \displaystyle Z_{\rm V}\;\bigg{[}\Big{(}1+3\overline{b}_{\rm V}\;am_{\rm q}^{% \rm av}+\frac{b_{\rm V}}{3}\;a(m_{{\rm q},l}+2m_{{\rm q},s})\Big{)}\,V_{\mu}^{% 8,I} | (4) | ||

\displaystyle\qquad\leavevmode\nobreak\ +\,({\textstyle\frac{1}{3}}b_{\rm V}+f% _{\rm V})\;\frac{2}{\sqrt{3}}a(m_{{\rm q},l}-m_{{\rm q},s})\;V_{\mu}^{0,I}% \bigg{]}\,, |

with V_{\mu}^{0}=\frac{1}{2}\bar{\psi}\gamma_{\mu}\psi the flavor-singlet current. Here V_{\mu}^{3,I} and V_{\mu}^{8,I} are understood to already contain the improvement term proportional to c_{\rm V}. All coefficients appearing in the two equations above are functions of the coupling \widetilde{g}_{0}. In this article, we present a non-perturbative determination of the renormalization factors Z_{\rm V}, b_{\rm V} and \overline{b}_{\rm V} as well as of c_{\rm V}, while the coefficient f_{\rm V} will remain undetermined. As explained below, there are however good reasons to expect f_{\rm V} to be numerically very small Bhattacharya:2005rb (). The improvement coefficient c_{\rm V} is determined by imposing continuum chiral Ward identities, as proposed in quenched QCD in Ref. Guagnelli:1997db (). We follow the presentation of Ref. Bhattacharya:2005rb () for the full renormalization and improvement in large volumes with N_{\mathrm{f}}=2+1 Wilson fermions. The mass-dependent renormalization with N_{\mathrm{f}}=2 Wilson fermions has been computed in Ref. Bakeyev:2003ff (). Note that the method of Ref. Harris:2015vfa () 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 \pi^{+}\pi^{-} 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 Z_{\rm V}, b_{\rm V} and \overline{b}_{\rm V}.

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_{\rm A} of the axial current, and appendix B contains some results on the employed correlation functions in lattice perturbation theory.

## II Renormalization and improvement: theory background

### II.1 Definitions and notations

We use Euclidean Dirac matrices, \{\gamma_{\mu},\gamma_{\nu}\}=2\delta_{\mu\nu}. We consider initially the general case of N_{\rm f} flavours of quarks. Flavor indices will be denoted by latin letters i,j,\dots Let

A_{\mu}^{(ij)}(x)=\overline{\psi}_{i}(x)\gamma_{\mu}\gamma_{5}\psi_{j}(x)\,,% \quad P^{(ij)}(x)=\overline{\psi}_{i}(x)\gamma_{5}\psi_{j}(x)\, | (5) |

be the bare axial current and pseudoscalar density. The on-shell improved operators are given by

(A_{I}^{(ij)})_{\mu}(x)=A_{\mu}^{(ij)}(x)+ac_{\rm A}(g_{0}^{2})\,\partial_{\mu% }P^{(ij)}(x)\,,\quad P_{I}^{(ij)}(x)=P^{(ij)}(x)\qquad(i\neq j)\,, | (6) |

where c_{\rm A} is an improvement coefficient. The average bare PCAC quark mass m_{ij} of quark flavours i and j is defined through the relation

\langle\partial_{\mu}(A^{(ij)}_{I})_{\mu}(x)\,P^{(ji)}(y)\rangle=2m_{ij}\,% \langle P_{I}^{(ij)}(x)\,P^{(ji)}(y)\rangle+{\rm O}(a^{2})\qquad(i\neq j,x\neq y% )\,. | (7) |

We also defined the subtracted bare quark mass of flavour i,

m_{{\rm q},i}=m_{0,i}-m_{\rm cr}. | (8) |

Often, the hopping parameter \kappa_{i}\equiv(8+2am_{0,i})^{-1} is used to parametrize the bare quark mass m_{0,i}. The value \kappa_{\rm cr}\equiv(8+2am_{\rm cr})^{-1} of the hopping parameter is the value for which the PCAC mass, defined through Eq. (7), vanishes in the {\rm SU}(N_{\rm f})-symmetric theory. The bare quark mass matrix is defined as M_{0}={\rm diag}(m_{0,1},\dots,m_{0,N_{\rm f}}), and similarly for the subtracted bare quark mass matrix, M_{\rm q}={\rm diag}(m_{{\rm q},1},\dots,m_{{\rm q},N_{\rm f}}). Finally, we also introduce the average quark masses

m_{{\rm q},ij}={\textstyle\frac{1}{2}}(m_{{\rm q},i}+m_{{\rm q},j})\,,\quad m_% {\rm q}^{\rm av}=\frac{1}{N_{\rm f}}\sum_{i=1}^{N_{\rm f}}m_{\rm q,i}\,. | (9) |

Here we will be concerned with the improvement and renormalization of the vector current V_{\mu}^{(ij)} on the lattice. Two discretizations are in common use, the local (l) and the point-split (c) lattice vector currents,

\displaystyle V_{\mu}^{l,(ij)}(x) | \displaystyle=\overline{\psi}_{i}(x)\gamma_{\mu}\psi_{j}(x)\,, | (10a) | ||

\displaystyle V_{\mu}^{c,(ij)}(x) | \displaystyle=\frac{1}{2}\left(\overline{\psi}_{i}(x+a\hat{\mu})(1+\gamma_{\mu% })U^{{\dagger}}_{\mu}(x)\psi_{j}(x)-\overline{\psi}_{i}(x)(1-\gamma_{\mu})U_{% \mu}(x)\psi_{j}(x+a\hat{\mu})\right)\,. | (10b) |

Instead of the point-split vector current, we actually consider the symmetrized version (cs)

V_{\mu}^{cs,(ij)}(x)=\frac{1}{2}\left(V_{\mu}^{c,(ij)}(x)+V_{\mu}^{c,(ij)}(x-a% \hat{\mu})\right)\,, | (11) |

which has the same properties under space-time reflections as the local vector current^{1}^{1}1The authors would like to thank Stefan Sint for pointing out this fact Frezzotti:2001ea (). This ensures that the same counterterms are present to remove O(a) artifacts

(V_{I}^{(ij)})_{\mu}(x)=V_{\mu}^{(ij)}(x)+ac_{\rm V}(g_{0}^{2})\,\widetilde{% \nabla}_{\nu}\Sigma_{\mu\nu}^{(ij)}(x)\,, | (12) |

with the local tensor current defined as

\Sigma_{\mu\nu}^{(ij)}=-\frac{1}{2}\overline{\psi}_{i}[\gamma_{\mu},\gamma_{% \nu}]\psi_{j}\,, | (13) |

and where we use the symmetric lattice derivative,

\widetilde{\nabla}_{\nu}\phi(x)=\frac{\phi(x+a\hat{\nu})-\phi(x-a\hat{\nu})}{2% a}\,. | (14) |

Generically, the renormalization pattern of the quark bilinears, including O(a) mass-dependent effects, has been derived in Ref. Bhattacharya:2005rb (). For the vector current, and writing V_{\mu} as a flavor matrix, it reads

\displaystyle\mathrm{tr}(\lambda V_{\mu})_{R} | \displaystyle= | \displaystyle Z_{\rm V}(\tilde{g}_{0}^{2})\Big{[}\left(1+N_{\rm f}\,\overline{% b}_{\rm V}(g_{0}^{2})\,am_{\rm q}^{\rm av}\right)\mathrm{tr}(\lambda V_{\mu}^{% I})+\frac{1}{2}b_{\rm V}(g_{0}^{2})\,\mathrm{tr}(\{\lambda,aM_{\rm q}\}V_{\mu}% ^{I}) | (15) | ||

\displaystyle\qquad\qquad+\;f_{\rm V}(g_{0}^{2})\,\mathrm{tr}(\lambda\,aM_{\rm q% })\,\mathrm{tr}(V_{\mu}^{I})\Big{]}\,, |

where

\tilde{g}_{0}^{2}\equiv g_{0}^{2}(1+b_{\rm g}\,am_{\rm q}^{\rm av}) | (16) |

is the modified bare coupling, which is in one-to-one correspondence with the lattice spacing, irrespective of the quark masses Luscher:1996sc (). The symbol ‘\mathrm{tr}’ refers to the trace over flavor indices and the \lambda is any element of the SU(N_{\rm f}) Lie algebra. The improvement coefficients c_{\rm V}, b_{\rm V}, \overline{b}_{\rm V} and f_{\rm V} are function of the bare coupling only; Z_{\rm V} has no anomalous dimension and does not depend on the renormalization scale.

Given that the coefficient b_{\rm g} is so far only known perturbatively, it is worth noting the following. If one Taylor-expands the function Z_{\rm V} and only keeps terms up to O(a), the expression (15) is equivalent to replacing the argument of Z_{\rm V} by g_{0}^{2} and then substituting \overline{b}_{\rm V} by

\overline{b}_{\rm V}^{\,\rm eff}(g_{0}^{2})\equiv\overline{b}_{\rm V}(g_{0}^{2% })+{\frac{1}{N_{\rm f}}}b_{\rm g}(g_{0}^{2})\,\frac{g_{0}^{2}}{Z_{\rm V}}\frac% {dZ_{\rm V}}{dg_{0}^{2}}. | (17) |

Therefore, the renormalization conditions we use for the vector current are only able to determine the combination \overline{b}_{\rm V}^{\,\rm eff}. In a second step, using the perturbative estimate of b_{\rm g}, we obtain a value for \overline{b}_{\rm V}. In the future, when a non-perturbative determination of b_{\rm g} becomes available, the value of \overline{b}_{\rm V} can be updated.

In section II.2, we describe the strategy used to determine the renormalization constant Z_{\rm V} and the improvement coefficients b_{\rm V}, \overline{b}_{\rm V}^{\,\rm eff} and c_{\rm V}.

### II.2 Vector Ward identities and determination of Z_{\rm V}, b_{\rm V} and \overline{b}_{\rm V}

We define an infinitesimal local vector transformation by

\displaystyle\delta\psi(y) | \displaystyle=\lambda\,\alpha(y)\,\psi(y)\,, | (18a) | ||

\displaystyle\delta\overline{\psi}(y) | \displaystyle=-\overline{\psi}(y)\alpha(y)\,\lambda\,, | (18b) |

where the matrix \lambda 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,

\Big{\langle}\frac{\delta\mathcal{O}}{\delta\alpha(y)}\Big{\rangle}=\Big{% \langle}\mathcal{O}\frac{\delta S}{\delta\alpha(y)}\Big{\rangle}\,, | (19) |

where S is the Euclidean action and \mathcal{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 Bochicchio:1985xa ()

\Big{\langle}\frac{\delta\mathcal{O}}{\delta\alpha(y)}\Big{\rangle}=a^{4}{% \nabla}_{\mu}^{*}\Big{\langle}{\rm tr}\{\lambda^{\intercal}\,V_{\mu}^{c}(y)\}% \,\mathcal{O}\Big{\rangle}+a^{4}\Big{\langle}\overline{\psi}(y)\left[M_{0},\,% \lambda\right]\psi(y)\,\mathcal{O}\Big{\rangle}\,, | (20) |

where {\nabla}_{\mu}^{*}\phi(y)=(\phi(y)-\phi(y-a\hat{\mu}))/a is the backward lattice derivative in the \mu-direction, \lambda^{\intercal} denotes the matrix transpose of \lambda and V_{\mu}^{c,(ij)}(y) corresponds to the point-split vector current defined in Eq. (10b). Using an operator \mathcal{O} with support which does not contain the site y and for \left[M,{\lambda}\right]=0, one simply recovers the conservation equation for the point-split vector current.

Working in components, we now consider the vector transformation

\displaystyle\delta\psi_{i}(y)=+\alpha(y)\psi_{i}(y)\,,\quad\delta\overline{% \psi}_{i}(y)=-\alpha(y)\overline{\psi}_{i}(y) | (21) |

for one specific flavor i. Then, using \mathcal{O}(x,z)=P^{(ji)}(x)P^{(ij)}(z) as a probe operator with i\neq j, one finds

\frac{\delta\mathcal{O}(x,z)}{\delta\alpha(y)}=\mathcal{O}(x,z)\delta(y-x)-% \mathcal{O}(x,z)\delta(y-z)\,, | (22) |

such that Eq. (20) reads

\Big{\langle}P^{(ji)}(x)P^{(ij)}(z)\Big{\rangle}\left(\delta(y-z)-\delta(y-x)% \right)=a^{4}{\nabla}_{y,\mu}^{*}\Big{\langle}P^{(ji)}(x)V_{\mu}^{c,(ii)}(y)P^% {(ij)}(z)\Big{\rangle}\,. | (23) |

Summing over the spatial vector \boldsymbol{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 \langle a^{3}\sum_{\boldsymbol{y}}P^{(ji)}(x)V_{0}^{c,(ii)}(y)P^{(ij)}(z)\rangle, 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

R(x_{0}-z_{0},y_{0}-z_{0})=\frac{\langle a^{6}\sum_{\boldsymbol{x},\boldsymbol% {y}}P^{(ji)}(x)V_{0}^{c,(ii)}(y)P^{(ij)}(z)\rangle}{\langle a^{3}\sum_{% \boldsymbol{x}}P^{(ji)}(x)P^{(ij)}(z)\rangle}\,, | (24) |

is unity such that the point-split vector current does not need any renormalization factor : Z_{\rm V}^{c}=1 and \overline{b}_{\rm V}^{c}=b_{\rm V}^{c}=f_{\rm V}^{c}=0. On the other hand, the local vector current is not conserved on the lattice and needs to be renormalized.

In N_{\rm f}=2+1 QCD with a quark mass matrix given by (2), by imposing either of the ratios

\displaystyle R_{\pi}(x_{0}-z_{0},y_{0}-z_{0}) | \displaystyle=\frac{\langle a^{6}\sum_{\boldsymbol{x},\boldsymbol{y}}P^{(21)}(% x)\,\frac{1}{2}(V_{0,\mathrm{R}}^{l,(11)}(y)-V_{0,\mathrm{R}}^{l,(22)}(y))P^{(% 12)}(z)\rangle}{\langle a^{3}\sum_{\boldsymbol{x}}P^{(21)}(x)P^{(12)}(z)% \rangle}\,, | (25a) | ||

\displaystyle R_{K}(x_{0}-z_{0},y_{0}-z_{0}) | \displaystyle=\frac{\langle a^{6}\sum_{\boldsymbol{x},\boldsymbol{y}}P^{(31)}(% x)\,(V_{0,\mathrm{R}}^{l,(11)}(y)-V_{0,\mathrm{R}}^{l,(22)}(y))P^{(13)}(z)% \rangle}{\langle a^{3}\sum_{\boldsymbol{x}}P^{(31)}(x)P^{(13)}(z)\rangle} | (25b) |

to equal unity on the lattice at finite quark masses, one can determine the renormalization factor of the local isovector current V_{\mu}^{3}=\frac{1}{2}(V_{\mu}^{(11)}-V_{\mu}^{(22)}), 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_{\rm V}, and that the two choices for the “spectator quark” correspond to two different renormalization prescriptions. Using ensembles with different values of m_{\rm q,1}=m_{{\rm q},2} and m_{\rm q,3}=m_{{\rm q},s}, each parameter can be determined independently. We remark that Z_{\rm V}, b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff} could also be determined in the same way from the matrix element of

\widetilde{R}_{K}(x_{0}-z_{0},y_{0}-z_{0})=\frac{-\langle a^{6}\sum_{% \boldsymbol{x},\boldsymbol{y}}P^{(31)}(x)\,\frac{1}{3}(V_{0,\mathrm{R}}^{l,(11% )}(y)+V_{0,\mathrm{R}}^{l,(22)}(y)-2V_{0,\mathrm{R}}^{l,(33)}(y))P^{(13)}(z)% \rangle}{\langle a^{3}\sum_{\boldsymbol{x}}P^{(31)}(x)P^{(13)}(z)\rangle}, | (26) |

since the flavor-singlet charge operator does not contribute on a kaon state. On the other hand, to obtain sensitivity to the coefficient f_{\rm V}, an external state with non-vanishing baryon number is required, for instance one may require the ratio

R_{\Delta^{++}}(x_{0}-z_{0},y_{0}-z_{0})=\frac{\langle a^{6}\sum_{\boldsymbol{% x},\boldsymbol{y}}\Delta^{(111)}(x)\,(V_{0,\mathrm{R}}^{l,(22)}(y)-V_{0,% \mathrm{R}}^{l,(33)}(y))\bar{\Delta}^{(111)}(z)\rangle}{\langle a^{3}\sum_{% \boldsymbol{x}}\Delta^{(111)}(x)\bar{\Delta}^{(111)}(z)\rangle} | (27) |

to vanish.
Without the vector current improvement term proportional to f_{\rm V}, R_{\Delta^{++}} would receive a contribution
of order a from disconnected diagrams; the role of the coefficient f_{\rm V}, which multiplies the flavor-singlet vector current, under which the \Delta^{++} baryon is charged, is to cancel this contribution. Therefore, the magnitude of f_{\rm V} is determined by the size of disconnected diagrams with the insertion of a single vector current. In perturbation theory, f_{\rm V} is therefore of order g_{0}^{6}, because at least three gluons must be emitted from the quark loop^{2}^{2}2One-gluon exchange does not contribute because the color factor vanishes.
To see that the two-gluon exchange also vanishes, one may use the \gamma_{5}-hermiticity of the quark propagator, \gamma_{5}S(x,y)\gamma_{5}=S(y,x)^{\dagger}, the fact that the free quark propagator S(x,y) is Hermitian for fixed (x,y) and \gamma_{5}\gamma_{\mu}\gamma_{5}=-\gamma_{\mu}, to show that the two orientations with which the quark propagators contribute to the quark loop come with opposite sign and cancel each other.

### II.3 Axial Ward identities and determination of c_{\rm V}

Once the renormalization factor Z_{\rm V} and improvement coefficients b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff} are known, the improvement coefficient c_{\rm V} can be determined by enforcing an axial Ward identity. In the continuum theory, the latter can be derived from the specific transformation

\delta\psi_{1}(x)=-\alpha(x)\gamma_{5}\psi_{2}(x)\,,\qquad\delta\overline{\psi% }_{2}(x)=-\overline{\psi}_{1}(x)\alpha(x)\gamma_{5}\,. | (28) |

As the operator to be chirally rotated, we choose {\cal O}(y)=A_{\mu}^{(23)}(y), and we have

\delta A_{\mu}^{(23)}(y)=\alpha(y)V_{\mu}^{(13)}(y)\,. | (29) |

Choosing \alpha(x) to be unity inside the slab t_{1}<x_{0}<t_{2} and zero outside, the variation of the action (per unit \alpha) is given by

\delta S^{(12)}=-\int_{t_{1}}^{t_{2}}\mathrm{d}x^{0}\,\int\mathrm{d}^{3}x\left% (2m_{R,12}P_{R}^{(12)}(x)-\partial_{\mu}A_{R,\mu}^{(12)}(x)\right)\,, | (30) |

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}\notin[t_{1},t_{2}], we then obtain the following Ward identity

\int\mathrm{d}^{3}y\,\Big{\langle}\delta S^{(12)}A_{R,k}^{(23)}(y_{0},% \boldsymbol{y})\,\mathcal{O}_{\rm ext}^{(31)}(z_{0},\boldsymbol{0})\Big{% \rangle}=\int\mathrm{d}^{3}y\,\Big{\langle}V_{R,k}^{(13)}(y_{0},\boldsymbol{y}% )\,\mathcal{O}_{\rm ext}^{(31)}(z_{0},\boldsymbol{0})\Big{\rangle}\,, | (31) |

valid in the continuum Guagnelli:1997db (), 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 (\alpha=l,cs), we define our estimator

\hat{c}^{\alpha}_{\mathrm{V}}(m_{q},g_{0}^{2})=\frac{\langle\sum_{\boldsymbol{% y}}\delta S^{(12)}A_{R,k}^{(23)}(y_{0},\boldsymbol{y})\,\mathcal{O}_{\rm ext}^% {(31)}(z_{0},\boldsymbol{0})-\hat{Z}_{V}^{(13)}\,\sum_{\boldsymbol{y}}V_{k}^{% \alpha,(13)}(y_{0},\boldsymbol{y})\,\mathcal{O}_{\rm ext}^{(31)}(z_{0},% \boldsymbol{0})\rangle}{\hat{Z}_{V}^{(13)}\langle\sum_{\boldsymbol{y}}a\,% \partial_{\nu}\Sigma_{k\nu}^{(13)}(y_{0},\boldsymbol{y})\,\mathcal{O}_{\rm ext% }^{(31)}(z_{0},\boldsymbol{0})\rangle} | (32) |

with

\hat{Z}_{V}^{(13)}(g_{0}^{2},m_{\rm q}^{\rm av},m_{\rm q,13})=Z_{\rm V}(g_{0}^% {2})\left(1+3\overline{b}_{\rm V}^{\,\rm eff}(g_{0}^{2})am_{\rm q}^{\rm av}+b_% {\rm V}(g_{0}^{2})am_{\rm q,13}\right) | (33) |

for the local vector current and \hat{Z}_{V}^{(13)}(g_{0}^{2},m_{\rm q}^{\rm av},m_{\rm 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_{\rm V}. The renormalized axial and pseudoscalar operators for i\neq j are respectively given by

\displaystyle P_{R}^{(ij)}(x) | \displaystyle=Z_{\rm P}\,(1+3\overline{b}_{\rm P}\,am_{\rm q}^{\rm av}+b_{\rm P% }\,am_{{\rm q},ij})\ P_{I}^{(ij)}(x)\,, | (34a) | ||

\displaystyle A_{k,R}^{(ij)}(x) | \displaystyle=Z_{\rm A}\,(1+3\overline{b}_{\rm A}\,am_{\rm q}^{\rm av}+b_{\rm A% }\,am_{{\rm q},ij})\,A_{k,I}^{(ij)}(x)\,, | (34b) |

in terms of the improved operators defined in Eq. (6). The renormalized quark mass is defined through the relation

m_{ij}=m_{R,ij}\,\frac{Z_{\rm P}(1+3\overline{b}_{P}\,am_{\rm q}^{\rm av}+b_{% \rm P}\,am_{{\rm q},ij})}{Z_{\rm A}(1+3\overline{b}_{\rm A}\,am_{\rm q}^{\rm av% }+b_{\rm A}\,am_{{\rm q},ij})}\qquad(i\neq j)\,, | (35) |

such that the combination m_{R,12}P_{R}^{(12)} is insensitive to Z_{\rm P}, b_{\rm P} and \overline{b}_{\rm P}. The renormalization factor Z_{\rm A} and the improvement parameters b_{\rm A} and \overline{b}_{\rm A} have been determined non-perturbatively in Refs. DallaBrida:2018tpn (); Bulava:2016ktf (); Korcyl:2016ugy ().

In Eq. (32), the operator \mathcal{O}_{\rm ext} can be either the vector \mathcal{O}_{\rm ext}=V_{k}^{(31)}(z_{0},\boldsymbol{0}) or the tensor operator \mathcal{O}_{\rm ext}=\Sigma_{k0}^{(31)}(z_{0},\boldsymbol{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}\to 0 must be taken to remove this contact term. This is done by computing the effective \hat{c}_{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_{\rm A} using an axial Ward identity. Our non-perturbative determination of c_{\rm A}, which we can compare to the literature Bulava:2015bxa (), serves as a cross-check of our numerical setup.

### II.4 Known perturbative results

The known perturbative results in QCD with N_{\rm c} colors and N_{\rm f} flavors of quarks are the following. The result b_{\rm g}=0.012000(2)\,N_{\rm f}\,g_{0}^{2}+{\rm O}(g_{0}^{4}), independently of the pure gauge action, was obtained in Sint:1995ch (). For degenerate quarks, only the combination b_{\rm V}+N_{\rm f}\overline{b}_{\rm V} appears, and the perturbative series for \overline{b}_{\rm V} starts at order g_{0}^{4}. For the tree-level improved Lüscher-Weisz action, the results are (C_{\rm F}=\frac{N_{\rm c}^{2}-1}{2N_{\rm c}}) Aoki:1998qd (); Taniguchi:1998pf ()

\displaystyle Z_{\rm V} | \displaystyle= | \displaystyle 1-0.075427\times C_{\rm F}\,g_{0}^{2}+{\rm O}(g_{0}^{4})\,, | (36a) | ||

\displaystyle b_{\rm V} | \displaystyle= | \displaystyle 0.0886(2)\times C_{\rm F}\,g_{0}^{2}+{\rm O}(g_{0}^{4})\,, | (36b) | ||

\displaystyle c_{\rm V}^{l} | \displaystyle= | \displaystyle-0.01030(4)\times C_{\rm F}\,g_{0}^{2}+{\rm O}(g_{0}^{4})\,. | (36c) |

The tree-level value of c_{\rm V}^{cs} is \frac{1}{2}.

Trajectory | CLS | \quad\beta\quad | L^{3}\times T | \kappa_{l} | \kappa_{s} | m_{\pi}\leavevmode\nobreak\ [\mathrm{MeV}] | am_{\rm PCAC} | \hat{Z}_{V}^{(12)} | \hat{Z}_{V}^{(12)} | \hat{c}_{V}^{l} | \hat{c}_{V}^{cs} | \hat{c}_{A} |

m_{q}^{\rm av}= const. with OBC | H101 | 3.40 | 32^{3}\times 96 | 0.13675962 | 0.13675962 | 420 | 0.009196(47) | 0.71562(4) | 0.71562(4) | +0.041(32) | +0.468(23) | -0.0404(25) |

H102 | 32^{3}\times 96 | 0.136865 | 0.13654934 | 350 | 0.006512(49) | 0.71226(6) | 0.71252(5) | -0.018(30) | +0.425(22) | -0.0505(20) | ||

H105 | 32^{3}\times 96 | 0.136970 | 0.13634079 | 280 | 0.004021(55) | 0.70908(5) | 0.70956(4) | -0.017(30) | +0.426(22) | -0.0508(24) | ||

N101 | 48^{3}\times 128 | 0.136970 | 0.13634079 | 290 | 0.003994(32) | 0.70911(5) | 0.72802(5) | \times | \times | \times | ||

C101 | 48^{3}\times 96 | 0.137030 | 0.13622204 | 220 | 0.002494(41) | 0.70717(9) | 0.70776(4) | \times | \times | \times | ||

S400 | 3.46 | 32^{3}\times 128 | 0.136984 | 0.13670239 | 350 | 0.005673(36) | 0.72355(5) | 0.72372(5) | +0.029(18) | +0.464(13) | -0.0354(26) | |

N401 | 48^{3}\times 128 | 0.1370616 | 0.13654808 | 290 | 0.003781(26) | 0.72116(5) | 0.72148(5) | -0.004(15) | +0.443(13) | -0.0406(21) | ||

H200 | 3.55 | 32^{3}\times 96 | 0.137000 | 0.137000 | 420 | 0.006858(23) | 0.74028(5) | 0.74028(5) | +0.021(30) | +0.458(22) | -0.0279(22) | |

N202 | 48^{3}\times 128 | 0.137000 | 0.137000 | 410 | 0.006856(16) | 0.74028(4) | 0.74028(4) | \times | \times | \times | ||

N203 | 48^{3}\times 128 | 0.137080 | 0.13684028 | 340 | 0.004761(13) | 0.73792(4) | 0.73802(5) | -0.018(29) | +0.431(23) | -0.0327(29) | ||

N200 | 48^{3}\times 128 | 0.137140 | 0.13672086 | 280 | 0.003145(15) | 0.73614(6) | 0.73632(4) | -0.002(28) | +0.442(22) | -0.0303(25) | ||

D200 | 64^{3}\times 128 | 0.137200 | 0.13660175 | 200 | 0.001538(12) | 0.73429(4) | 0.73461(4) | \times | \times | \times | ||

N300 | 3.70 | 48^{3}\times 128 | 0.137000 | 0.137000 | 420 | 0.005509(09) | 0.75909(2) | 0.75909(2) | -0.047(25) | +0.409(20) | -0.0205(18) | |

N302 | 48^{3}\times 128 | 0.137064 | 0.13687218 | 340 | 0.003712(11) | 0.75719(3) | 0.75728(3) | -0.063(26) | +0.399(20) | -0.0213(25) | ||

J303 | 64^{3}\times 192 | 0.137123 | 0.13675466 | 260 | 0.002047(08) | 0.75544(2) | 0.75557(3) | -0.036(28) | +0.418(24) | -0.0206(33) | ||

m_{q,s}=m_{q,s}^{\rm phys} with OBC | H107 | 3.40 | 32^{3}\times 96 | 0.13694567 | 0.13620317 | 360 | 0.006731(81) | 0.71056(8) | 0.71110(6) | \times | \times | \times |

H106 | 32^{3}\times 96 | 0.13701557 | 0.13614870 | 280 | 0.004111(73) | 0.70792(12) | 0.70853(6) | \times | \times | \times | ||

C102 | 48^{3}\times 96 | 0.13705085 | 0.13612906 | 230 | 0.002573(42) | 0.70676(8) | 0.70731(6) | \times | \times | \times | ||

N204 | 3.55 | 48^{3}\times 128 | 0.137112 | 0.13657505 | 350 | 0.004884(26) | 0.73715(4) | 0.73745(4) | \times | \times | \times | |

N201 | 48^{3}\times 128 | 0.13715968 | 0.13656132 | 280 | 0.003166(34) | 0.73564(6) | 0.73600(6) | \times | \times | \times | ||

D201 | 64^{3}\times 128 | 0.1372067 | 0.13654684 | 200 | 0.001526(25) | 0.73401(5) | 0.73436(5) | \times | \times | \times | ||

N304 | 3.70 | 48^{3}\times 128 | 0.13707933 | 0.13666543 | 350 | 0.003603(43) | 0.75683(6) | 0.75702(6) | \times | \times | \times | |

m_{q,l}=m_{q,s} with PBC | rqcd29 | 3.46 | 32^{3}\times 64 | 0.136600 | 0.136600 | 710 | 0.021687(77) | 0.73741(6) | 0.73741(6) | \times | \times | \times |

B450 | 32^{3}\times 64 | 0.136890 | 0.136890 | 410 | 0.008071(36) | 0.72647(2) | 0.72647(2) | +0.020(16) | +0.455(12) | -0.0403(11) | ||

rqcd30 | 32^{3}\times 64 | 0.1369587 | 0.1369587 | 320 | 0.004832(59) | 0.72398(6) | 0.72398(6) | \times | \times | \times | ||

X450 | 48^{3}\times 64 | 0.136994 | 0.136994 | 260 | 0.003305(31) | 0.72265(3) | 0.72265(3) | \times | \times | \times | ||

X250 | 3.55 | 48^{3}\times 64 | 0.137050 | 0.137050 | 350 | 0.004919(23) | 0.73863(2) | 0.73863(2) | \times | \times | \times | |

X251 | 48^{3}\times 64 | 0.137100 | 0.137100 | 270 | 0.002877(26) | 0.73698(2) | 0.73698(2) | \times | \times | \times | ||

m_{q,l}=m_{q,s} with OBC | H401 | 3.46 | 32^{3}\times 96 | 0.136725 | 0.136725 | 590 | 0.015694(67) | 0.73260(8) | 0.73260(8) | \times | \times | \times |

H400 | 32^{3}\times 96 | 0.13688848 | 0.13688848 | 420 | 0.008235(81) | 0.72634(10) | 0.72634(10) | \times | \times | \times | ||

N303 | 3.70 | 48^{3}\times 128 | 0.136800 | 0.136800 | 650 | 0.012518(20) | 0.76544(2) | 0.76544(2) | \times | \times | \times | |

N301 | 48^{3}\times 128 | 0.137005 | 0.137005 | 420 | 0.005334(19) | 0.75897(3) | 0.75897(3) | \times | \times | \times |

## III Numerical setup

\beta | T/(2a) | z_{0}/a | [t_{1},t_{2}]/a | y_{0}/a |
---|---|---|---|---|

3.40 | 48 | 41 | [46,54] | 49-50 |

3.46 | 32 | 0 | [6,15] | 10-11 |

48 | 0 | [6,15] | 10-11 | |

3.55 | 48 | 41 | [48,59] | 52-53 |

64 | 57 | [64,75] | 68-69 | |

3.70 | 64 | 53 | [62,76] | 68-69 |

We use the N_{f}=2+1 CLS (Coordinated Lattice Simulations) lattice ensembles whose main parameters are given in Table. 1.
They have been produced using the openQCD code^{3}^{3}3http://luscher.web.cern.ch/luscher/openQCD/ of Ref. Luscher:2012av () using the Wilson-Clover action for fermions and the tree-level Symanzik-improved gauge action. The parameter c_{\rm SW} has been determined non-perturbatively in Ref. Bulava:2013cta (). We consider four values of the bare coupling \beta=3.40,3.46,3.55 and 3.70 which correspond to lattice spacings in the range 0.050-0.085\leavevmode\nobreak\ \mathrm{fm} Bruno:2016plf (). 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_{\rm q}^{\rm av} and m_{{\rm q},s}=m_{{\rm q},s}^{\rm phys} can be used 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_{{\rm q},l}=m_{{\rm q},s}.
Concerning c_{\rm V}, it is enough to consider ensembles on a single chiral trajectory (e.g. m_{\rm q}^{\rm av}= constant). However, to determine the
two improvement coefficients b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff}, we have to consider at least two different chiral trajectories.

For the calculation of the renormalization factor Z_{\rm V}, we need to compute the following three-point correlation function, projected on vanishing momentum

C_{PVP}(x_{0},y_{0};z_{0})=a^{6}\sum_{\boldsymbol{x},\boldsymbol{y}}\langle P^% {(ij)}(x_{0},\boldsymbol{x})V^{(jj)}_{0}(y_{0},\boldsymbol{y})P^{(ij){\dagger}% }(z_{0},\boldsymbol{0})\rangle\,, | (37) |

and two-point correlation functions

\displaystyle C_{PP}(x_{0},z_{0}) | \displaystyle= | \displaystyle a^{3}\sum_{\boldsymbol{x}}\langle P^{(ij)}(x_{0},\boldsymbol{x})% P^{(ij){\dagger}}(z_{0},\boldsymbol{0})\rangle\,, | (38a) | ||

\displaystyle C_{AP}(x_{0},z_{0}) | \displaystyle= | \displaystyle a^{3}\sum_{\boldsymbol{x}}\langle A^{(ij)}_{0}(x_{0},\boldsymbol% {x})P^{(ij){\dagger}}(z_{0},\boldsymbol{0})\rangle\,. | (38b) |

Correlation functions are calculated using U(1) stochastic sources with time dilution Foley:2005ac (). On each gauge configuration, we generate an ensemble of N_{s} stochastic sources with support on a single timeslice and satisfying

\lim_{N_{s}\rightarrow\infty}\frac{1}{N_{s}}\sum_{s=1}^{N_{s}}\eta_{\alpha}^{a% }(x)_{s}\left[\eta_{\beta}^{b}(y)_{s}\right]^{*}=a^{-3}\delta_{\alpha\beta}% \delta^{ab}\delta_{x,y}\,, | (39) |

where each component is normalized to one, \eta_{\alpha}^{a}(x)_{[r]}^{*}\ \eta_{\alpha}^{a}(x)_{[r]}=1 (no summation). This can be implemented by using U(1) noise for each color and spinor index on site x of the lattice. For ensembles with open boundary conditions in the time direction, time translation is lost. In this case, the source is placed at z_{0}=T/4 away from the boundary (t=0) and the two-point correlation functions are obtained for all values of x_{0}\in[0,T/2], keeping the sink time away from the second boundary (t=T). For the three-point correlation function, the sink time is placed at x_{0}=3T/4 and is computed for all y_{0}.

For each stochastic source s with support in time-slice z_{0}, we solve the Dirac equation and denote the solution vector \Phi_{i}^{s}(x;z_{0})=a^{3}\sum_{\boldsymbol{z}}S(x,z)\eta_{i}^{s}(z). Correlation functions are given by

\displaystyle C^{(ij)}_{PP}(x_{0},z_{0}) | \displaystyle=\frac{a^{6}}{V}\sum_{\boldsymbol{x},\boldsymbol{z}}{\rm Tr}\left% [S_{i}(x,z)^{{\dagger}}S_{j}(x,z)\right]=\frac{a^{3}}{N_{s}V}\sum_{s,% \boldsymbol{x}}\Phi_{i}^{s}(x;z_{0})^{{\dagger}}\Phi^{s}_{j}(x;z_{0})\,, | (40a) | ||

\displaystyle C^{(ij)}_{AP}(x_{0},z_{0}) | \displaystyle=\frac{a^{6}}{V}\sum_{\boldsymbol{x},\boldsymbol{z}}{\rm Tr}\left% [S_{i}(x,z)^{{\dagger}}\gamma_{0}S_{j}(x,z)\right]=\frac{a^{3}}{N_{s}V}\sum_{% \boldsymbol{x}}\Phi^{s}_{i}(x;z_{0})^{{\dagger}}\gamma_{0}\Phi^{s}_{j}(x;z_{0}% )\,, | (40b) | ||

\displaystyle C^{(ij)}_{PVP}(x_{0},y_{0};z_{0}) | \displaystyle=\frac{a^{9}}{V}\sum_{\boldsymbol{x},\boldsymbol{y},\boldsymbol{z% }}{\rm Tr}\left[S_{i}(x,z)^{{\dagger}}S_{j}(x,y)\gamma_{0}S_{j}(y,z)\right] | (40c) | ||

\displaystyle=\frac{a^{3}}{N_{s}V}\sum_{\boldsymbol{x}}\widetilde{\Phi}_{ji}^{% s{\dagger}}(y;x_{0},z_{0})\gamma_{0}\Phi^{s}_{j}(x;z_{0})\,, |

with V=L^{3} the spatial volume. We have used the \gamma_{5}-hermiticity of the fermion propagator S(x,y)=\gamma_{5}S(y,x)^{{\dagger}}\gamma_{5} and \widetilde{\Phi}^{s}_{ji} is a sequential propagator given by \widetilde{\Phi}^{s}_{ji}(y;x_{0},z_{0})=a^{6}\sum_{\boldsymbol{x},\boldsymbol%
{z}}\gamma_{5}S_{j}(y,x)\gamma_{5}S_{i}(x,z)\eta_{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 2. 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

### IV.1 Results for Z_{\rm V}, b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff}

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 \hat{Z}_{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

R(t,t_{1})\xrightarrow[t_{1},\;t-t_{1}\to\infty]{}\frac{1}{\hat{Z}_{V}^{(12)}}\,, | (41) |

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_{PVP}(x_{0},y_{0};z_{0})\to C_{PVP}(x_{0},y_{0};z_{0})-C_{PVP}(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 1. In a second step, the renormalization constant Z_{\rm V}, and the improvement coefficients b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff} at a given value of the bare coupling g_{0} are obtained using the fit ansatz

\hat{Z}_{\rm V}^{(12)}(g_{0}^{2},m_{\rm q}^{\rm av},m_{\rm q,12})=Z_{\rm V}(g_% {0}^{2})\left(1+3\overline{b}_{\rm V}^{\,\rm eff}(g_{0}^{2})am_{\rm q}^{\rm av% }+b_{\rm V}(g_{0}^{2})am_{\rm q,12}\right)\,. | (42) |

The ensembles included in the fit satisfy |am_{q,l}|<0.01 and am_{\rm q}^{\rm av}<0.01 such that higher order corrections are expected to be small. The results for each value of \beta are given in Table 3 and the statistical error includes the error on \kappa_{\rm cr}. We note that the coefficient b_{\rm V} is significantly larger than the one-loop perturbative estimate given in Eq. (36b) and that \overline{b}_{\rm V}^{\,\rm eff}\ll b_{\rm V}. This was expected since the perturbative value starts only at two loops in perturbation theory. We provide the covariance matrices for the different values of the coupling considered in this work

\displaystyle\mathrm{cov}(Z_{\rm V},b_{\rm V},\overline{b}_{\rm V}^{\,\rm eff}% ;\beta=3.40) | \displaystyle=\begin{pmatrix}+6.04\times 10^{-8}&-1.05\times 10^{-6}&-5.02% \times 10^{-6}\\ -1.05\times 10^{-6}&+1.93\times 10^{-4}&+1.30\times 10^{-4}\\ -5.02\times 10^{-6}&+1.30\times 10^{-4}&+5.50\times 10^{-4}\\ \end{pmatrix}\,, | (43a) | ||

\displaystyle\mathrm{cov}(Z_{\rm V},b_{\rm V},\overline{b}_{\rm V}^{\,\rm eff}% ;\beta=3.46) | \displaystyle=\begin{pmatrix}+4.16\times 10^{-9}&+9.92\times 10^{-8}&-6.68% \times 10^{-8}\\ +9.92\times 10^{-8}&+1.90\times 10^{-4}&-6.64\times 10^{-5}\\ -6.68\times 10^{-8}&-6.64\times 10^{-5}&+2.97\times 10^{-5}\\ \end{pmatrix}\,, | (43b) | ||

\displaystyle\mathrm{cov}(Z_{\rm V},b_{\rm V},\overline{b}_{\rm V}^{\,\rm eff}% ;\beta=3.55) | \displaystyle=\begin{pmatrix}+3.17\times 10^{-9}&-2.85\times 10^{-7}&-1.27% \times 10^{-7}\\ -2.85\times 10^{-7}&+9.21\times 10^{-5}&+1.37\times 10^{-5}\\ -1.27\times 10^{-7}&+1.37\times 10^{-5}&+1.37\times 10^{-5}\\ \end{pmatrix}\,, | (43c) | ||

\displaystyle\mathrm{cov}(Z_{\rm V},b_{\rm V},\overline{b}_{\rm V}^{\,\rm eff}% ;\beta=3.70) | \displaystyle=\begin{pmatrix}+3.60\times 10^{-9}&+2.10\times 10^{-7}&-1.48% \times 10^{-7}\\ +2.10\times 10^{-7}&+1.38\times 10^{-4}&-6.14\times 10^{-5}\\ -1.48\times 10^{-7}&-6.14\times 10^{-5}&+3.26\times 10^{-5}\\ \end{pmatrix}\,. | (43d) |

Finally, we perform a Padé fit to obtain the renormalization factor and the improvement coefficients as a function of the bare coupling g_{0}^{2}. Our final results read

\displaystyle Z_{\rm V}(g_{0}^{2}) | \displaystyle=1-0.10057\,g_{0}^{2}\times\frac{1+p_{1}\,g_{0}^{2}+p_{2}\,g_{0}^% {4}}{1+p_{3}\,g_{0}^{2}}\,, | (44a) | ||

\displaystyle b_{\rm V}(g_{0}^{2}) | \displaystyle=1+0.11813\,g_{0}^{2}\times\frac{1+p_{1}\,g_{0}^{2}}{1+p_{2}\,g_{% 0}^{2}}\,, | (44b) | ||

\displaystyle\overline{b}_{\rm V}^{\,\rm eff}(g_{0}^{2}) | \displaystyle=\frac{p_{1}\,g_{0}^{4}}{1+p_{2}\,g_{0}^{2}}\,, | (44c) |

which automatically reproduce the one-loop calculations and where the parameters and covariance matrices are given by

\displaystyle p(Z_{\rm V}) | \displaystyle=\begin{pmatrix}-0.2542\\ -0.0961\\ -0.4796\\ \end{pmatrix}\,,\qquad\mathrm{cov}(Z_{\rm V})=\begin{pmatrix}+1.31619&+4.92750% &+6.15758\\ +4.92750&+66.8321&+75.3218\\ +6.15758&+75.3218&+85.2733\\ \end{pmatrix}\times 10^{-6}\,, | (45a) | ||

\displaystyle p(b_{\rm V}) | \displaystyle=\begin{pmatrix}-0.184\\ -0.444\\ \end{pmatrix}\,,\qquad\ \ \ \mathrm{cov}(b_{\rm V})=\begin{pmatrix}+36.7139&+1% 2.6698\\ +12.6698&+4.41224\\ \end{pmatrix}\times 10^{-4}\,, | (45b) | ||

\displaystyle p(\overline{b}_{\rm V}^{\,\rm eff}) | \displaystyle=\begin{pmatrix}+0.00112\\ -0.5577\\ \end{pmatrix}\,,\qquad\mathrm{cov}(\overline{b}_{\rm V}^{\,\rm eff})=\begin{% pmatrix}+1.061463&+14.53004\\ +14.53004&+248.5266\\ \end{pmatrix}\times 10^{-8}\,. | (45c) |

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 \hat{Z}^{(12)}_{\rm 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 1 are in perfect agreement within statistical errors. Secondly, the correlation functions in Eq. (24) are computed with a source located at z_{0}=T/4 to suppress boundary effects. For the ensemble N101, we have performed three sets of simulations with different source locations, z_{0}=T/4, T/4-4a, T/4-8a, and the results are \hat{Z}^{(12)}_{\rm V}=0.70910(6), 0.70911(5) and 0.70919(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_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff} 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 3. We do not observe any difference for the renormalization constant Z_{\rm V} at our level of precision (both results should differ only by an O(a^{2}) ambiguity). For b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff}, 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_{\rm V}, we observe the expected linear scaling with the lattice spacing. For \overline{b}_{\rm V}, the ratio goes to one only if one includes higher order discretization effects, which appear to be sizeable.

#### IV.1.1 Comparison of results with previous work

The renormalization factor Z_{\rm V} has been determined independently in DallaBrida:2018tpn () using the chirally rotated Schrödinger functional framework. In Fig. 4, we plot the ratio Z_{\rm V,\,sub}^{l}/Z_{\rm V} where Z_{\rm V,\,sub}^{l} is extracted from DallaBrida:2018tpn () 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 \beta=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_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff} have also been determined recently in Ref. Fritzsch:2018zym () using a different setup, based on the QCD Schrödinger functional. A comparison with our results is given in Fig 6. For b_{\rm 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 \overline{b}_{\rm V}^{\,\rm eff}, the difference with the results quoted in Ref. Fritzsch:2018zym () is significant, especially at large couplings g_{0}^{2}. 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_{\rm q}^{\rm av} is 0.005 on the \beta=3.55 ensembles, so that with 3\overline{b}_{\rm V}^{\,\rm eff}\simeq 0.16 even a 100\% ambiguity on \overline{b}_{\rm 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 \overline{b}_{\rm V}^{\,\rm eff}. For the physics applications discussed in the introduction, it is interesting to compare our values for the renormalization factor \hat{Z}^{(12)}_{\rm V} of the isovector current to those one would obtain using the Z_{\rm V} values of DallaBrida:2018tpn () and the b_{\rm V} and \overline{b}_{\rm V}^{\,\rm eff} values from Fritzsch:2018zym (). We find that our direct estimates are always slightly lower, and that the relative difference depends almost only on g_{0}^{2}: it is about 1.3% at \beta=3.40, 0.8% at \beta=3.46, 0.37% at \beta=3.55 and 0.12% at \beta=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.

### IV.2 Results for the improvement coefficient c_{\rm V}

For y_{0} not too close from t_{1} and t_{2}, we can extract the value of \hat{c}_{\rm 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\leavevmode\nobreak\ \mathrm{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 2. 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 \mathcal{O}_{\rm ext}=V_{k}^{(31)}(z_{0},\boldsymbol{0})+\Sigma_{k0}^{(31)}(z_% {0},\boldsymbol{0}) which helps to improve the statistical precision. For Z_{\rm A}, we use the results called Z^{l}_{\rm A,sub} using the L_{1}-LCP from DallaBrida:2018tpn (). For b_{\rm A} we used the published values in Korcyl:2016ugy () and for \overline{b}_{\rm A} we use values from private () (in practice, we use \overline{b}_{\rm A}^{\rm eff} which includes the b_{\rm g}-term for Z_{\rm A}, as in Eq. (17)). The results for \hat{c}_{\rm V} for each ensemble are given in Table 1 and differ from c_{\rm V} by the presence of a contact term which vanishes in the chiral limit, as explained in Sec. II.3. The improvement coefficient c_{\rm V} is obtained using a linear extrapolation in the light-quark PCAC mass m_{12} at constant m_{\rm q}^{\rm av} and the results, for each value of the bare coupling, are summarized in Table 3. 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 \overline{b}_{\rm A}. The uncertainty on Z_{\rm A} and b_{\rm A} appears to have a negligible impact at our level of precision. Finally, we perform linear or quadratic fits in g_{0}^{2} to determine c_{\rm V} as a function of the bare coupling. The results for the local and the conserved vector currents read

\displaystyle c_{\rm V}^{l}(g_{0}^{2}) | \displaystyle=-0.01030\,C_{\rm F}\,g_{0}^{2}\,\times\,(1-0.15(35)\,g_{0}^{2})\,, | (46a) | ||

\displaystyle c_{\rm V}^{cs}(g_{0}^{2}) | \displaystyle=\frac{1}{2}\,\times\,(1-0.093(13)\,g_{0}^{2})\,. | (46b) |

The parametrization is consistent with the perturbative predictions collected in Section II.4. Our values for c_{\rm V}^{l} are significantly smaller (in magnitude) than the preliminary values determined in Heitger:2017njs () 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_{\rm V}. However, our values for Z_{\rm V} also differ by more that one standard deviation from the ones computed in Heitger:2017njs (). Since Z_{\rm V} enters in the determination of c_{\rm V}, this could partly explain the disagreement. For example, if we use the preliminary results for Z_{\rm V} given in Heitger:2017njs (), our value would be \hat{c}_{\rm V}^{l}=-0.146 for H105 and \hat{c}_{\rm V}^{l}=-0.178 for N200. These observations highlight the need for a good control over Z_{\rm V} to determine precisely c_{\rm V}^{l}. Since the point-split vector current is conserved, Z_{\rm V}^{c}=1, this issue is absent for the determination of c_{\rm V}^{cs}.

## 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_{\rm f}=2+1 O(a)-improved Wilson quark action and the tree-level Symanzik-improved gauge action. Only one coefficient, f_{\rm V}, is missing but is also expected to be small, as it starts at O(g_{0}^{6}) in perturbation theory; in this regard, we note that for the two other mass-dependent improvement coefficients, b_{\rm V} and \overline{b}_{\rm 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 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_{\rm V}(\widetilde{g}_{0}) at the value of the renormalized coupling \widetilde{g}_{0} instead of the bare coupling g_{0}. We have taken this difference into account by replacing the improvement coefficient \overline{b}_{\rm V} by the effective parameter \overline{b}_{\rm V}^{\,\rm eff}, thus avoiding the use of the unknown coefficient b_{\rm 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_{\rm V}, b_{\rm V} and \overline{b}_{\rm V}.

Our final results for the different \beta values used in CLS simulations are summarized in Table 3. We also provide interpolating formulae through Eqs. (44) and (46). As a cross-check of our methods, we have recomputed the improvement coefficient c_{\rm A} and find good agreement with the results of Ref. Bulava:2015bxa (), 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_{\rm V}^{\alpha} 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. Fritzsch:2018zym () shows that a potentially-large O(a)-ambiguity in \overline{b}_{\rm 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. DallaBrida:2018tpn (). However, we note the relative discrepancy is rather small, and smaller than that observed for two different normalization conditions for the axial current DallaBrida:2018tpn (); Bulava:2016ktf ().

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_{\rm f}=2+1 hadronic contribution to the running of the weak mixing angle involves the flavor-singlet vector current, whose improvement coefficient \bar{c}_{\rm V} is unknown. A method to determine the latter based on a chiral Ward identity was proposed in Bhattacharya:2005rb ().

\beta | 3.40 | 3.46 | 3.55 | 3.70 |
---|---|---|---|---|

\kappa_{\rm cr} | 0.1369115(27) | 0.1370645(10) | 0.1371726(13) | 0.1371576(8) |

Z_{\rm V} | 0.70908(25) | 0.71998(6) | 0.73454(6) | 0.75413(6) |

0.70912(19) | 0.71998(6) | 0.73453(6) | 0.75413(6) | |

b_{\rm V} | 1.648(14) | 1.622(14) | 1.541(10) | 1.488(12) |

1.546(10) | 1.526(13) | 1.460(09) | 1.427(12) | |

\overline{b}_{\rm V}^{\,\rm eff} | 0.206(23) | 0.108(05) | 0.053(04) | 0.029(06) |

0.240(17) | 0.140(05) | 0.081(04) | 0.049(06) | |

\overline{b}_{\rm V} | 0.227(23) | 0.125(05) | 0.067(04) | 0.040(06) |

0.260(17) | 0.157(05) | 0.095(04) | 0.060(06) | |

c_{\rm V}^{l} | -0.069(32) | -0.008(20) | -0.031(32) | -0.039(29) |

c_{\rm V}^{cs} | 0.389(23) | 0.438(16) | 0.422(26) | 0.416(24) |

###### Acknowledgements.

We thank Stefan Sint, Rainer Sommer and Mattia Dalla Brida for helpful discussions and Gunnar Bali and Piotr Korcyl for sharing with us preliminary results of the improvement coefficients of the axial current. We are grateful for the access to the ensembles used here, made available to us through CLS. Correlation functions were computed on the platforms “Clover” at the Helmholtz-Institut Mainz and “Mogon II” at Johannes Gutenberg University Mainz. This work was partly supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme through grant agreement no. 771971-SIMDAMA. The CLS lattice ensembles used here were partly generated through computing time provided by the Gauss Centre for Supercomputing (GCS) through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC).## Appendix A Determination of the axial improvement coefficient c_{\rm A}

In this appendix, we use a similar setup to determine the improvement coefficient c_{\rm A} of the axial vector current (see Eq. (6)). The latter was previously determined non-perturbatively in Bulava:2016ktf () in the framework of the Schrödinger functional. This study can be seen as a consistency check of our method.

Within our numerical setup, described in Sec III, we can replace the axial vector current and the external operator \mathcal{O}_{\mathrm{ext}} in Eq. (31) by any other operator without new inversion of the Dirac operator. Therefore, we consider the following axial Ward identity

\int\mathrm{d}^{3}y\,\langle\delta S^{(12)}V_{R,0}^{(23)}(y_{0},\boldsymbol{y}% )\,\mathcal{O}_{\rm ext}^{(31)}(z_{0},\boldsymbol{0})\rangle=\int\mathrm{d}^{3% }y\,\langle A_{R,0}^{(13)}(y_{0},\boldsymbol{y})\,\mathcal{O}_{\rm ext}^{(31)}% (z_{0},\boldsymbol{0})\rangle\,, | (47) |

with \mathcal{O}_{\rm ext}^{(31)}=P^{(31)}. The variation of the action is given by Eq. (30), and, similarly to Eq. (31), with the constraint y_{0}\notin[t_{1},t_{2}]. Here, V_{R,0}^{(23)} and A_{R,0}^{(13)} are the renormalized and O(a)-improved vector and axial vector currents defined in Eq. (6). If one knows Z_{\rm V}, b_{\rm V}, \overline{b}_{\rm V}^{\,\rm eff}, b_{\rm A} and \overline{b}_{\rm A}, then c_{\rm A} can be determined by imposing this equation to hold on the lattice up to O(a^{2}) discretization effects, as done in Eq. (32) for c_{\rm V}.

Using the same procedure as for c_{\rm V} and with the local vector current, we obtain the results summarized in Table 1 for a subset of the ensembles. Similarly to c_{\rm V}, the chiral dependence is very mild and the contribution from the contact term in the left hand side of Eq. (47) is removed by taking the limit m_{12}\to 0. We obtain the results quoted in Table 4. As can be seen on Fig. 7, our results are close to the values quoted in Bulava:2016ktf () obtained using a different method. Since this improvement coefficient has an O(a)-ambiguity, we can attribute this small difference to the different schemes used.

\beta | 3.40 | 3.46 | 3.55 | 3.70 |
---|---|---|---|---|

c_{\rm A} | -0.060(4) | -0.038(4) | -0.033(5) | -0.021(3) |

## Appendix B Axial Ward identities in the free theory

In the following, we give the tree-level expressions in lattice perturbation theory for the correlation functions involved in the chiral Ward identities, Eqs. (31) and (47). We have used these expressions to test the lattice QCD code implementing the Ward identities.

We provide a more general expression in that we allow for a general spatial momentum, on the other hand we restrict ourselves to the equal-mass case, m_{1}=m_{2}=m_{3}=m. At order g_{0}^{0} with a Wilson quark action, the correlation functions do not depend on c_{\rm sw}. We use the standard notation

\displaystyle\hat{p}_{\mu}=\frac{2}{a}\sin\frac{ap_{\mu}}{2},\qquad\qquad p^{{% }^{\!\!\!\circ}}_{\mu}=\frac{1}{a}\sin ap_{\mu}, | (48) |

as well as

\displaystyle A(\boldsymbol{p}) | \displaystyle= | \displaystyle 1+am+\frac{1}{2}a^{2}\hat{\boldsymbol{p}}^{2}, | (49) | ||

\displaystyle B(\boldsymbol{p}) | \displaystyle= | \displaystyle m^{2}+(1+am)\hat{\boldsymbol{p}}^{2}+\frac{1}{2}a^{2}\sum_{k<l}% \hat{p}_{k}^{2}\hat{p}_{l}^{2}, | (50) | ||

\displaystyle C(\boldsymbol{p}) | \displaystyle= | \displaystyle m+\frac{a}{2}\Big{(}{\hat{\boldsymbol{p}}^{2}}-\frac{B(% \boldsymbol{p})}{A(\boldsymbol{p})}\Big{)}, | (51) | ||

\displaystyle{\cal D}_{\boldsymbol{p}} | \displaystyle= | \displaystyle\sqrt{B(\boldsymbol{p})\,(4A(\boldsymbol{p})+a^{2}B(\boldsymbol{p% }))}=\frac{2}{a}A(\boldsymbol{p})\sinh(a\omega_{\boldsymbol{p}}). |

Let \pm i\omega_{\boldsymbol{p}} be the pole in p_{0} of the fermion propagator. We note the identities

\displaystyle\frac{4}{a^{2}}\sinh^{2}(a\omega_{\boldsymbol{p}}/2) | \displaystyle= | \displaystyle\frac{B(\boldsymbol{p})}{A(\boldsymbol{p})}, | (52) | ||

\displaystyle\frac{1}{a^{2}}\sinh^{2}(a\omega_{\boldsymbol{p}}) | \displaystyle= | \displaystyle\boldsymbol{p^{{}^{\!\!\!\circ}}}{}^{2}+C(\boldsymbol{p})^{2}. | (53) |

The free fermion propagator in the time-momentum representation reads, with the convention {\rm sign}(0)=0,

\displaystyle S(x,y) | \displaystyle\equiv | \displaystyle\langle\psi(x)\bar{\psi}(y)\rangle=\int_{B}\frac{d^{3}\boldsymbol% {p}}{(2\pi)^{3}}\frac{e^{-\omega_{\boldsymbol{p}}|x_{0}-y_{0}|+i\boldsymbol{p}% \cdot(\boldsymbol{x}-\boldsymbol{y})}}{{\cal D}_{\boldsymbol{p}}} | ||

\displaystyle\left({\rm sgn}(x_{0}-y_{0})\frac{1}{a}\sinh(a\omega_{\boldsymbol% {p}})\gamma_{0}-i\boldsymbol{\gamma}\cdot\boldsymbol{p^{{}^{\!\!\!\circ}}}+C(% \boldsymbol{p})+\delta_{x_{0},y_{0}}\frac{1}{a}\sinh(a\omega_{\boldsymbol{p}})% \right), |

where \int_{B} denotes integration over the Brillouin zone, -\frac{\pi}{a}<p_{i}<\frac{\pi}{a}.

In all three-point functions in this appendix, we assume z_{0}<{\rm min}(y_{0},x_{0}). The correlation functions relevant to the Ward identity (47) are

\displaystyle a^{3}\sum_{\boldsymbol{y}}\;e^{-i\boldsymbol{p}\cdot(\boldsymbol% {y}-\boldsymbol{z})}\,\Big{\langle}A_{0}^{13}(y)\,P^{31}(z)\Big{\rangle} | (55) | ||

\displaystyle=4\,{\rm sign}(y_{0}-z_{0})\int_{B}\frac{d^{3}q}{(2\pi)^{3}}\,% \frac{C(\boldsymbol{q})\sinh\omega_{\boldsymbol{p}+\boldsymbol{q}}+C(% \boldsymbol{p}+\boldsymbol{q})\sinh(\omega_{\boldsymbol{q}})}{{\cal D}_{% \boldsymbol{q}}{\cal D}_{\boldsymbol{p}+\boldsymbol{q}}}e^{-|y_{0}-z_{0}|(% \omega_{\boldsymbol{p}+\boldsymbol{q}}+\omega_{\boldsymbol{q}})}, | |||

\displaystyle a^{6}\sum_{\boldsymbol{x},\boldsymbol{z}}\;e^{i\boldsymbol{p}% \cdot(\boldsymbol{z}-\boldsymbol{y})}\Big{\langle}A_{0}^{12}(x)\,V_{0}^{23}(y)% \,P^{31}(z)\Big{\rangle} | (56) | ||

\displaystyle=8\int_{B}\frac{d^{3}q}{(2\pi)^{3}}\,\frac{e^{-2\omega_{% \boldsymbol{q}}\theta(x_{0}-y_{0})(x_{0}-y_{0})+(z_{0}-y_{0})(\omega_{% \boldsymbol{p}+\boldsymbol{q}}+\omega_{\boldsymbol{q}})}}{{\cal D}_{% \boldsymbol{q}}^{2}{\cal D}_{\boldsymbol{p}+\boldsymbol{q}}}f(\boldsymbol{p},% \boldsymbol{q}), |

with

f(\boldsymbol{p},\boldsymbol{q})=\left\{\begin{array}[]{l@{~~~}l}C(\boldsymbol% {q}){\boldsymbol{q}}^{{}^{\!\!\!\circ}}\cdot{\boldsymbol{k}}^{{}^{\!\!\!\circ}% }-C(\boldsymbol{p}+\boldsymbol{q}){\boldsymbol{q}}^{{}^{\!\!\!\circ}}{}^{2}% \Big{|}_{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}\leavevmode\nobreak\ % \leavevmode\nobreak\ \leavevmode\nobreak\ &x_{0}<y_{0},\\ C(\boldsymbol{q})\Big{(}\sinh(\omega_{\boldsymbol{q}})\sinh(\omega_{% \boldsymbol{p}+\boldsymbol{q}})+{\boldsymbol{q}}^{{}^{\!\!\!\circ}}\cdot{% \boldsymbol{k}}^{{}^{\!\!\!\circ}}+C(\boldsymbol{q})C(\boldsymbol{p}+% \boldsymbol{q})\Big{)}_{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}% \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ &x_{0}>y_{0}% \end{array}\right. | (57) |

In Eq. (56), we have assumed x_{0}\neq y_{0}. Finally, under the same assumptions we obtain

\displaystyle G_{0}(x_{0}-z_{0},y_{0}-z_{0},\boldsymbol{p})\equiv a^{6}\sum_{% \boldsymbol{x},\boldsymbol{z}}e^{i\boldsymbol{p}\cdot(\boldsymbol{z}-% \boldsymbol{y})}\langle P^{12}(x)\;V^{23}_{0}(y)P^{31}(z)\rangle | (58) | ||

\displaystyle=-8\theta(x_{0}-y_{0})\int_{B}\frac{d^{3}q}{(2\pi)^{3}}\,\frac{e^% {-2\omega_{\boldsymbol{q}}\theta(x_{0}-y_{0})(x_{0}-y_{0})+(z_{0}-y_{0})(% \omega_{\boldsymbol{p}+\boldsymbol{q}}+\omega_{\boldsymbol{q}})}}{{\cal D}_{% \boldsymbol{q}}^{2}{\cal D}_{\boldsymbol{p}+\boldsymbol{q}}}\,g(\boldsymbol{p}% ,\boldsymbol{q}), |

with

g(\boldsymbol{p},\boldsymbol{q})=\sinh(\omega_{\boldsymbol{q}})\Big{(}\sinh(% \omega_{\boldsymbol{q}})\sinh(\omega_{\boldsymbol{k}})+{\boldsymbol{q}}^{{}^{% \!\!\!\circ}}\cdot{\boldsymbol{k}}^{{}^{\!\!\!\circ}}+C(\boldsymbol{q})C(% \boldsymbol{k})\Big{)}_{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}. | (59) |

The special case x_{0}=y_{0} must be treated separately,

\displaystyle G_{0}(y_{0}-z_{0},y_{0}-z_{0},\boldsymbol{p}) | \displaystyle= | \displaystyle-8\int_{B}\frac{d^{3}q}{(2\pi)^{3}}\frac{e^{-(y_{0}-z_{0})(\omega% _{\boldsymbol{p}+\boldsymbol{q}}+\omega_{\boldsymbol{q}})}}{{\cal D}_{% \boldsymbol{q}}^{2}{\cal D}_{\boldsymbol{p}+\boldsymbol{q}}} | ||

\displaystyle\times\frac{\sinh(\omega_{\boldsymbol{q}})}{2}\Big{(}[\sinh(% \omega_{\boldsymbol{q}})+C(\boldsymbol{q})][\sinh(\omega_{\boldsymbol{k}})+C(% \boldsymbol{k})]+{\boldsymbol{q}}^{{}^{\!\!\!\circ}}\cdot{\boldsymbol{k}}^{{}^% {\!\!\!\circ}}\Big{)}_{\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q}}. |

The correlation functions relevant to the Ward identity Eq. (31) are

\displaystyle a^{3}\sum_{\boldsymbol{y}}e^{i\boldsymbol{p}\cdot(\boldsymbol{z}% -\boldsymbol{y})}\Big{\langle}V_{3}^{l,13}(y)\;\Sigma_{30}^{31}(z)\Big{\rangle% }=-4\;{\rm sign}(y_{0}-z_{0}) | (61) | ||

\displaystyle\qquad\times\int_{B}\frac{d^{3}\boldsymbol{q}}{(2\pi)^{3}}\,\frac% {e^{-(\omega_{\boldsymbol{p}}+\omega_{\boldsymbol{p}+\boldsymbol{q}})|z_{0}-y_% {0}|}}{{\cal D}_{\boldsymbol{q}}{\cal D}_{\boldsymbol{p}+\boldsymbol{q}}}\,% \Big{(}C(\boldsymbol{q})\sinh\omega_{\boldsymbol{p}+\boldsymbol{q}}+C(% \boldsymbol{p}+\boldsymbol{q})\sinh\omega_{\boldsymbol{q}}\Big{)}, |

and, for x_{0}\neq y_{0} and \boldsymbol{k}_{\perp}=(k_{1},k_{2}), \boldsymbol{\ell}_{\perp}=(\ell_{1},\ell_{2}),

\displaystyle a^{6}\sum_{\boldsymbol{x},\boldsymbol{z}}e^{i\boldsymbol{p}\cdot% (\boldsymbol{z}-\boldsymbol{y})}\Big{\langle}A_{0}^{12}(x)\;A_{3}^{23}(y)\;% \Sigma_{30}^{31}(z)\Big{\rangle} | (62) | ||

\displaystyle=-16\int_{B}\frac{d^{3}k}{(2\pi)^{3}}\,\frac{e^{-2\omega_{% \boldsymbol{k}}(x_{0}-y_{0})\theta(x_{0}-y_{0})-(\omega_{\boldsymbol{k}}+% \omega_{\boldsymbol{p}-\boldsymbol{k}})(y_{0}-z_{0})}}{{\cal D}_{\boldsymbol{k% }}^{2}\,{\cal D}_{\boldsymbol{p}-\boldsymbol{k}}}\,f^{A}(\boldsymbol{p},% \boldsymbol{k}), | |||

\displaystyle f^{A}(\boldsymbol{p},\boldsymbol{k})=\frac{1}{2}\Big{\{}\sinh(% \omega_{\boldsymbol{k}})\Big{(}C(\boldsymbol{k})\sinh(\omega_{\boldsymbol{\ell% }})\theta(x_{0}-y_{0})-\sinh(\omega_{\boldsymbol{k}})C(\boldsymbol{\ell})% \theta(y_{0}-x_{0})\Big{)} | |||

\displaystyle-C(\boldsymbol{k})k^{{}^{\!\!\!\circ}}_{3}\ell^{{}^{\!\!\!\circ}}% _{3}+C(\boldsymbol{k}){\boldsymbol{k}}^{{}^{\!\!\!\circ}}_{\perp}\cdot{% \boldsymbol{\ell}}^{{}^{\!\!\!\circ}}_{\perp}+C(\boldsymbol{k})^{2}C(% \boldsymbol{\ell})\Big{\}}_{\boldsymbol{\ell}=\boldsymbol{p}-\boldsymbol{k}}, | (63) | ||

\displaystyle f^{A}(\boldsymbol{0},\boldsymbol{k})=\left\{\begin{array}[]{l@{~% ~~}l}C(\boldsymbol{k})(k^{{}^{\!\!\!\circ}}_{3}{}^{2}+C(\boldsymbol{k})^{2})% \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ &x_{0}>y_{0}\\ -C(\boldsymbol{k}){\boldsymbol{k}}^{{}^{\!\!\!\circ}}_{\perp}{}^{2}\leavevmode% \nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ &x_{0}<y_{0}\end{array}% \right.. | (64) |

Further, for x_{0}\neq y_{0}, we have

\displaystyle G^{A}(x_{0}-z_{0},y_{0}-z_{0},\boldsymbol{p})=a^{6}\sum_{% \boldsymbol{x},\boldsymbol{z}}e^{i\boldsymbol{p}\cdot(\boldsymbol{z}-% \boldsymbol{y})}\Big{\langle}P^{12}(x)A_{3}^{23}(y)\Sigma_{30}^{31}(z)\Big{\rangle} | (65) | ||

\displaystyle=-16\,\theta(x_{0}-y_{0})\int_{B}\frac{d^{3}k}{(2\pi)^{3}}\,\frac% {e^{-2\omega_{\boldsymbol{k}}(x_{0}-y_{0})\theta(x_{0}-y_{0})+(\omega_{% \boldsymbol{k}}+\omega_{\boldsymbol{p}-\boldsymbol{k}})(z_{0}-y_{0})}}{{\cal D% }_{\boldsymbol{k}}^{2}\,{\cal D}_{\boldsymbol{p}-\boldsymbol{k}}}\,g^{A}(% \boldsymbol{p},\boldsymbol{k}), | |||

\displaystyle g^{A}(\boldsymbol{p},\boldsymbol{k})=-\frac{1}{2}\sinh(\omega_{% \boldsymbol{k}})\Big{(}\sinh(\omega_{\boldsymbol{k}})\sinh(\omega_{\boldsymbol% {\ell}})-k^{{}^{\!\!\!\circ}}_{3}\ell^{{}^{\!\!\!\circ}}_{3}+{\boldsymbol{k}}^% {{}^{\!\!\!\circ}}_{\perp}\cdot{\boldsymbol{\ell}}^{{}^{\!\!\!\circ}}_{\perp}+% C(\boldsymbol{k})C(\boldsymbol{\ell})\Big{)}_{\boldsymbol{\ell}=\boldsymbol{p}% -\boldsymbol{k}}, | (66) | ||

\displaystyle g^{A}(\boldsymbol{0},\boldsymbol{k})=-\sinh(\omega_{\boldsymbol{% k}})(k^{{}^{\!\!\!\circ}}_{3}{}^{2}+C(\boldsymbol{k})^{2}). | (67) |

Finally, for x_{0}=y_{0},

\displaystyle G^{A}(y_{0}-z_{0},y_{0}-z_{0},\boldsymbol{p})=-4\int_{B}\frac{d^% {3}k}{(2\pi)^{3}}\,\frac{e^{(\omega_{\boldsymbol{k}}+\omega_{\boldsymbol{p}-% \boldsymbol{k}})(z_{0}-y_{0})}}{{\cal D}_{\boldsymbol{k}}^{2}\,{\cal D}_{% \boldsymbol{p}-\boldsymbol{k}}}\,\sinh\omega_{\boldsymbol{k}} | (68) | ||

\displaystyle\times\Big{[}k^{{}^{\!\!\!\circ}}_{3}\ell^{{}^{\!\!\!\circ}}_{3}-% {\boldsymbol{k}}^{{}^{\!\!\!\circ}}_{\perp}\cdot{\boldsymbol{\ell}}^{{}^{\!\!% \!\circ}}_{\perp}-C_{\boldsymbol{\ell}}C_{\boldsymbol{k}}-\sinh\omega_{% \boldsymbol{k}}\sinh\omega_{\boldsymbol{\ell}}-C_{\boldsymbol{\ell}}\sinh% \omega_{\boldsymbol{k}}-C_{\boldsymbol{k}}\sinh\omega_{\boldsymbol{\ell}}\Big{% ]}_{\boldsymbol{\ell}=\boldsymbol{p}-\boldsymbol{k}}. |

## References

- (1) S. Aoki et al., Eur. Phys. J. C 77, no. 2, 112 (2017) [arXiv:1607.00299 [hep-lat]].
- (2) H. B. Meyer and H. Wittig, arXiv:1807.09370 [hep-lat].
- (3) T. Bhattacharya, R. Gupta, W. Lee, S. R. Sharpe and J. M. S. Wu, Phys. Rev. D 73, 034504 (2006) [hep-lat/0511014].
- (4) M. Guagnelli and R. Sommer, Nucl. Phys. Proc. Suppl. 63, 886 (1998) [hep-lat/9709088].
- (5) T. Bakeyev et al. [QCDSF-UKQCD Collaboration], Phys. Lett. B 580, 197 (2004) [hep-lat/0305014].
- (6) T. Harris and H. B. Meyer, Phys. Rev. D 92, 114503 (2015) [arXiv:1506.05248 [hep-lat]].
- (7) R. Frezzotti et al. [ALPHA Collaboration], JHEP 0107, 048 (2001) [hep-lat/0104014].
- (8) M. Luscher, S. Sint, R. Sommer and P. Weisz, Nucl. Phys. B 478, 365 (1996) [hep-lat/9605038].
- (9) M. Bochicchio, L. Maiani, G. Martinelli, G. C. Rossi and M. Testa, Nucl. Phys. B 262, 331 (1985).
- (10) P. Korcyl and G. S. Bali, Phys. Rev. D 95, no. 1, 014505 (2017) [arXiv:1607.07090 [hep-lat]].
- (11) J. Bulava, M. Della Morte, J. Heitger and C. Wittemeier, Phys. Rev. D 93, no. 11, 114513 (2016) [arXiv:1604.05827 [hep-lat]].
- (12) M. Dalla Brida, T. Korzec, S. Sint and P. Vilaseca, arXiv:1808.09236 [hep-lat].
- (13) S. Sint and R. Sommer, Nucl. Phys. B 465 (1996) 71 [hep-lat/9508012].
- (14) S. Aoki, R. Frezzotti and P. Weisz, Nucl. Phys. B 540 (1999) 501 [hep-lat/9808007].
- (15) Y. Taniguchi and A. Ukawa, Phys. Rev. D 58 (1998) 114503 [hep-lat/9806015].
- (16) M. Bruno, T. Korzec and S. Schaefer, Phys. Rev. D 95, no. 7, 074504 (2017) [arXiv:1608.08900 [hep-lat]].
- (17) M. Lüscher and S. Schaefer, Comput. Phys. Commun. 184, 519 (2013) [arXiv:1206.2809 [hep-lat]].
- (18) J. Bulava and S. Schaefer, Nucl. Phys. B 874, 188 (2013) [arXiv:1304.7093 [hep-lat]].
- (19) J. Foley, K. Jimmy Juge, A. O’Cais, M. Peardon, S. M. Ryan and J. I. Skullerud, Comput. Phys. Commun. 172, 145 (2005) [hep-lat/0505023].
- (20) P. Fritzsch, JHEP 1806, 015 (2018) [arXiv:1805.07401 [hep-lat]].
- (21) P. Korcyl and G. S. Bali, private communication.
- (22) S. Sint and P. Weisz, Nucl. Phys. B 502, 251 (1997) [hep-lat/9704001].
- (23) J. Heitger, F. Joswig, A. Vladikas and C. Wittemeier, arXiv:1711.03924 [hep-lat].
- (24) J. Bulava et al. [ALPHA Collaboration], Nucl. Phys. B 896, 555 (2015) [arXiv:1502.04999 [hep-lat]].
- (25) G. S. Bali et al. [RQCD Collaboration], Phys. Rev. D 94 (2016) no.7, 074501 [arXiv:1606.09039 [hep-lat]].