Non-perturbative renormalization and O (a) -improvement of the non-singlet vector current with N f =2+1 Wilson fermions and tree-level Symanzik improved gauge action
MMITP/18-106
Non-perturbative renormalization and O( a )-improvement of the non-singletvector current with N f = 2 + 1 Wilson fermions and tree-level Symanzikimproved gauge action Antoine G´erardin,
1, 2
Tim Harris,
3, 4 and Harvey B. Meyer
1, 3 Institut f¨ur Kernphysik & Cluster of Excellence PRISMA,Johannes Gutenberg-Universit¨at Mainz, D-55099 Mainz, Germany John von Neumann Institute for Computing,DESY, Platanenallee 6, D-15738 Zeuthen, Germany Helmholtz Institut Mainz, D-55099 Mainz, Germany Dipartimento di Fisica G. Occhialini, Universita’ di Milano-Bicocca, and,INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy (Dated: February 20, 2019)In calculating hadronic contributions to precision observables for tests of the StandardModel in lattice QCD, the electromagnetic current plays a central role. Using a Wilsonaction with O( a ) improvement in QCD with N f flavors, a counterterm must be added to thevector current in order for its on-shell matrix elements to be O( a ) improved. In addition, thelocal vector current, which has support on one lattice site, must be renormalized. At O( a ),the breaking of the SU( N f ) symmetry by the quark mass matrix leads to a mixing betweenthe local currents of different quark flavors. We present a non-perturbative calculation ofall the required improvement and renormalization constants needed for the local and theconserved electromagnetic current in QCD with N f = 2 + 1 O( a )-improved Wilson fermionsand tree-level Symanzik improved gauge action, with the exception of one coefficient, whichwe show to be order g in lattice perturbation theory. The method is based on the vectorand axial Ward identities imposed at finite lattice spacing and in the chiral limit. We makeuse of lattice ensembles generated as part of the Coordinated Lattice Simulations (CLS)initiative. I. INTRODUCTION
Precision tests of the Standard Model typically require reliable theory input from first-principles calculations. While the electroweak sector can be treated perturbatively, the virtualcontributions 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 ofthe Standard Model [1]. Example of such hadronic quantities are the ratio of decay constants f K /f π , the MS quark masses, the running strong coupling constant α s ( M Z ), and the anomalousmagnetic moment of the muon, ( g − µ . For the latter, a major effort by several lattice collab-orations worldwide is ongoing to calculate the hadronic vacuum polarization and the hadroniclight-by-light contributions [2]. From the QCD point of view, these contributions amount totwo- and four-point correlation functions of the electromagnetic current, to be integrated overwith a weight function containing the characteristic scale of the muon mass.In continuum QCD, the electromagnetic current is conserved and does not require renormal-ization. On the lattice, a finite renormalization can appear, depending on the details of theaction 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 thecurrent, which is exactly conserved at finite lattice spacing. This appealing property, however,does not guarantee that transverse correlation functions of the current have smaller discretiza-tion effects than those of the naive, “local” vector current with support on a single lattice site,which in the limit of massless quarks receives a finite renormalization factor Z V ( g ). Indeed,the improvement of the vector current – local or point-split – requires adding the divergence ofthe tensor current with a coefficient denoted c V , which counteracts the breaking of chiral sym-metry by the Wilson action and suffices to remove all O( a ) cutoff effects in on-shell correlationfunctions. This coefficient, whose value depends on the discretization of the current, has a finitevalue at tree level of perturbation theory in the case of the point-split current, but vanishes forthe local current.On the other hand, for the local vector current, a mass dependence of the renormaliza- a r X i v : . [ h e p - l a t ] F e b tion factor arises if O( a ) discretization errors are to be removed. This mass dependence isrelevant in precision calculations, given the pattern of the physical up, down and strangequark masses. Concretely, given the electric charge matrix of the lightest three quark flavors, Q = diag(2 / , − / , − / V e . m .µ = V µ + 1 √ V µ , (1)where V aµ = ¯ ψγ µ λ a ψ is the octet of vector currents, with λ a the Gell-Mann matrices. In isospin-symmetric QCD, the bare quark mass matrix can be decomposed as M q = m avq + 1 √ m q ,l − m q ,s ) λ . (2)See Eq. (8) and below for our notation. The renormalization pattern of the local discretizationof the two neutral octet combinations then reads [3], at O( a ), V µ,R = Z V (1 + 3 b V am avq + b V am q ,l ) V ,Iµ , (3) V µ,R = Z V (cid:20)(cid:16) b V am avq + b V a ( m q ,l + 2 m q ,s ) (cid:17) V ,Iµ + ( b V + f V ) 2 √ a ( m q ,l − m q ,s ) V ,Iµ (cid:21) , (4)with V µ = ¯ ψγ µ ψ the flavor-singlet current. Here V ,Iµ and V ,Iµ are understood to alreadycontain the improvement term proportional to c V . All coefficients appearing in the two equa-tions above are functions of the coupling (cid:101) g . In this article, we present a non-perturbativedetermination of the renormalization factors Z V , b V and b V as well as of c V , while the coeffi-cient f V will remain undetermined. As explained below, there are, however, good reasons toexpect f V to be numerically very small [3]. The improvement coefficient c V is determined byimposing continuum chiral Ward identities, as proposed in quenched QCD in Ref. [4]. We followthe presentation of Ref. [3] for the full renormalization and improvement in large volumes with N f = 2 + 1 Wilson fermions. The mass-dependent renormalization with N f = 2 Wilson fermionshas been computed in Ref. [5]. Note that the method of Ref. [6] allows only one to compute alinear combination of the improvement coefficients for the conserved and local currents, and itis 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 theelectromagnetic current with only O( a ) discretization effects. This will in particular allow for ashorter continuum extrapolation of the leading hadronic contribution to the anomalous magneticmoment of the muon, and therefore a more cost-effective set of lattice QCD simulations. Giventhat phenomenologically, the π + π − channel, which is described by the timelike electromagneticform 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 theelectric charge of the pion be unity at every lattice spacing. This is the main renormalizationcondition we will adopt to determine Z V , b V and b V .We begin by giving an overview of the required theory background, which allows us to defineour notation. We present the setup for our numerical calculation in Sec. III and the results inSec. IV. We finish with some concluding remarks in Sec. V. Appendix A presents a determinationof the improvement coefficient c A of the axial current, and Appendix B contains some resultson the employed correlation functions in lattice perturbation theory. II. RENORMALIZATION AND IMPROVEMENT: THEORY BACKGROUNDA. Definitions and notations
We use Euclidean Dirac matrices, { γ µ , γ ν } = 2 δ µν . We consider initially the general case of N f flavors of quarks. Flavor indices will be denoted by latin letters i, j, . . . Let A ( ij ) µ ( x ) = ψ i ( x ) γ µ γ ψ j ( x ) , P ( ij ) ( x ) = ψ i ( x ) γ ψ j ( x ) (5)be the bare axial current and pseudoscalar density. The on-shell improved operators are givenby ( A ( ij ) I ) µ ( x ) = A ( ij ) µ ( x ) + ac A ( g ) ∂ µ P ( ij ) ( x ) , P ( ij ) I ( x ) = P ( ij ) ( x ) ( i (cid:54) = j ) , (6)where c A is an improvement coefficient. The average bare PCAC quark mass m ij of quark flavors i and j is defined through the relation (cid:104) ∂ µ ( A ( ij ) I ) µ ( x ) P ( ji ) ( y ) (cid:105) = 2 m ij (cid:104) P ( ij ) I ( x ) P ( ji ) ( y ) (cid:105) + O( a ) ( i (cid:54) = j, x (cid:54) = y ) . (7)We also defined the subtracted bare quark mass of flavor i , m q ,i = m ,i − m cr . (8)Often, the hopping parameter κ i ≡ (8 + 2 am ,i ) − is used to parametrize the bare quark mass m ,i . The value κ cr ≡ (8 + 2 am cr ) − of the hopping parameter is the value for which the PCACmass, defined through Eq. (7), vanishes in the SU( N f )-symmetric theory. The bare quark massmatrix is defined as M = diag( m , , . . . , m ,N f ), and similarly for the subtracted bare quarkmass matrix, M q = diag( m q , , . . . , m q ,N f ). Finally, we also introduce the average quark masses m q ,ij = ( m q ,i + m q ,j ) , m avq = 1 N f N f (cid:88) i =1 m q , i . (9)Here we will be concerned with the improvement and renormalization of the vector current V ( ij ) µ on the lattice. Two discretizations are in common use, the local ( l ) and the point-split ( c )lattice vector currents, V l, ( ij ) µ ( x ) = ψ i ( x ) γ µ ψ j ( x ) , (10a) V c, ( ij ) µ ( x ) = 12 (cid:16) ψ i ( x + a ˆ µ )(1 + γ µ ) U † µ ( x ) ψ j ( x ) − ψ i ( x )(1 − γ µ ) U µ ( x ) ψ j ( x + a ˆ µ ) (cid:17) . (10b)Instead of the point-split vector current, we actually consider the symmetrized version ( cs ) V cs, ( ij ) µ ( x ) = 12 (cid:16) V c, ( ij ) µ ( x ) + V c, ( ij ) µ ( x − a ˆ µ ) (cid:17) , (11)which has the same properties under spacetime reflections as the local vector current [7]. Thisensures that the same counterterms are present to remove O( a ) artifacts( V ( ij ) I ) µ ( x ) = V ( ij ) µ ( x ) + ac V ( g ) (cid:101) ∇ ν Σ ( ij ) µν ( x ) , (12)with the local tensor current defined asΣ ( ij ) µν = − ψ i [ γ µ , γ ν ] ψ j , (13)and where we use the symmetric lattice derivative, (cid:101) ∇ ν φ ( x ) = φ ( x + a ˆ ν ) − φ ( x − a ˆ ν )2 a . (14)Generically, the renormalization pattern of the quark bilinears, including O( a ) mass-dependenteffects, has been derived in Ref. [3]. For the vector current, and writing V µ as a flavor matrix,it reads tr( λV µ ) R = Z V (˜ g ) (cid:104) (cid:0) N f b V ( g ) am avq (cid:1) tr( λV Iµ ) + 12 b V ( g ) tr( { λ, aM q } V Iµ )+ f V ( g ) tr( λ aM q ) tr( V Iµ ) (cid:105) , (15) The authors thank Stefan Sint for pointing out this fact where ˜ g ≡ g (1 + b g am avq ) (16)is the modified bare coupling, which is in one-to-one correspondence with the lattice spacing,irrespective of the quark masses [8]. The symbol ‘tr’ refers to the trace over flavor indices and λ is any element of the SU( N f ) Lie algebra. The improvement coefficients c V , b V , b V and f V arefunctions of the bare coupling only; Z V has no anomalous dimension and does not depend onthe renormalization scale.Given that the coefficient b g is so far only known perturbatively, it is worth noting the follow-ing. If one Taylor expands the function Z V and only keeps terms up to O( a ), the expression (15)is equivalent to replacing the argument of Z V by g and then substituting b V by b effV ( g ) ≡ b V ( g ) + 1 N f b g ( g ) g Z V dZ V dg . (17)Therefore, the renormalization conditions we use for the vector current are only able to determinethe combination b effV . In a second step, using the perturbative estimate of b g , we obtain a valuefor b V . In the future, when a non-perturbative determination of b g becomes available, the valueof b V can be updated.In Sec. II B, we describe the strategy used to determine the renormalization constant Z V andthe improvement coefficients b V , b effV and c V . B. Vector Ward identities and determination of Z V , b V and b V We define an infinitesimal local vector transformation by δψ ( y ) = λ α ( y ) ψ ( y ) , (18a) δψ ( y ) = − ψ ( y ) α ( y ) λ , (18b)where the matrix λ acts on flavor space. Using the path integral definition of an expectationvalue and noticing that the previous transformation is a change of integration variables withunit Jacobian, one obtains the following identity, (cid:68) δ O δα ( y ) (cid:69) = (cid:68) O δSδα ( y ) (cid:69) , (19)where S is the Euclidean action and O is any operator. In fact, the equality holds on everysingle gauge-field configuration because only the fermionic part of the action is affected. ForWilson-Clover fermions, it leads to the well-known vector Ward identity [9] (cid:68) δ O δα ( y ) (cid:69) = a ∇ ∗ µ (cid:68) tr { λ (cid:124) V cµ ( y ) } O (cid:69) + a (cid:68) ψ ( y ) [ M , λ ] ψ ( y ) O (cid:69) , (20)where ∇ ∗ µ φ ( y ) = ( φ ( y ) − φ ( y − a ˆ µ )) /a is the backward lattice derivative in the µ -direction, λ (cid:124) denotes the matrix transpose of λ and V c, ( ij ) µ ( y ) corresponds to the point-split vector currentdefined in Eq. (10b). Using an operator O with support which does not contain the site y andfor [ M, λ ] = 0, one simply recovers the conservation equation for the point-split vector current.Working in components, we now consider the vector transformation δψ i ( y ) = + α ( y ) ψ i ( y ) , δψ i ( y ) = − α ( y ) ψ i ( y ) (21)for one specific flavor i . Then, using O ( x, z ) = P ( ji ) ( x ) P ( ij ) ( z ) as a probe operator with i (cid:54) = j ,one finds δ O ( x, z ) δα ( y ) = O ( x, z ) δ ( y − x ) − O ( x, z ) δ ( y − z ) , (22)such that Eq. (20) reads (cid:68) P ( ji ) ( x ) P ( ij ) ( z ) (cid:69) ( δ ( y − z ) − δ ( y − x )) = a ∇ ∗ y,µ (cid:68) P ( ji ) ( x ) V c, ( ii ) µ ( y ) P ( ij ) ( z ) (cid:69) . (23)Summing over the spatial vector y in Eq. (23), the spatial derivative does not contribute dueto the use of periodic boundary conditions and only the time derivative remains. Therefore,the three-point correlation function (cid:104) a (cid:80) y P ( ji ) ( x ) V c, ( ii )0 ( y ) P ( ij ) ( z ) (cid:105) , viewed as a function of y , is a piecewise constant function with discrete steps of +1 at y = z and − y = x . Inparticular, for x > y > z , the ratio R defined by R ( x − z , y − z ) = (cid:104) a (cid:80) x , y P ( ji ) ( x ) V c, ( ii )0 ( y ) P ( ij ) ( z ) (cid:105)(cid:104) a (cid:80) x P ( ji ) ( x ) P ( ij ) ( z ) (cid:105) , (24)is unity such that the point-split vector current does not need any renormalization factor : Z c V = 1 and b c V = b c V = f c V = 0. On the other hand, the local vector current is not conserved onthe lattice and needs to be renormalized.In N f = 2 + 1 QCD with a quark mass matrix given by (2), by imposing either of the ratios R π ( x − z , y − z ) = (cid:104) a (cid:80) x , y P (21) ( x ) ( V l, (11)0 , R ( y ) − V l, (22)0 , R ( y )) P (12) ( z ) (cid:105)(cid:104) a (cid:80) x P (21) ( x ) P (12) ( z ) (cid:105) , (25a) R K ( x − z , y − z ) = (cid:104) a (cid:80) x , y P (31) ( x ) ( V l, (11)0 , R ( y ) − V l, (22)0 , R ( y )) P (13) ( z ) (cid:105)(cid:104) a (cid:80) x P (31) ( x ) P (13) ( z ) (cid:105) (25b)to equal unity on the lattice at finite quark masses, one can determine the renormalization factorof the local isovector current V µ = ( V (11) µ − V (22) µ ), including the O( a ) mass-dependent terms,as given explicitly in Eq. (3). We note that this renormalization condition does not requirethe knowledge of c V , and that the two choices for the “spectator quark” correspond to twodifferent renormalization prescriptions. Using ensembles with different values of m q , = m q , and m q , = m q ,s , each parameter can be determined independently. We remark that Z V , b V and b effV could also be determined in the same way from the matrix element of (cid:101) R K ( x − z , y − z ) = −(cid:104) a (cid:80) x , y P (31) ( x ) ( V l, (11)0 , R ( y ) + V l, (22)0 , R ( y ) − V l, (33)0 , R ( y )) P (13) ( z ) (cid:105)(cid:104) a (cid:80) x P (31) ( x ) P (13) ( z ) (cid:105) , (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 V , an external state with nonvanishing baryon numberis required, for instance one may require the ratio R ∆ ++ ( x − z , y − z ) = (cid:104) a (cid:80) x , y ∆ (111) ( x ) ( V l, (22)0 , R ( y ) − V l, (33)0 , R ( y )) ¯∆ (111) ( z ) (cid:105)(cid:104) a (cid:80) x ∆ (111) ( x ) ¯∆ (111) ( z ) (cid:105) (27)to vanish. Without the vector current improvement term proportional to f V , R ∆ ++ wouldreceive a contribution of order a from disconnected diagrams; the role of the coefficient f V ,which multiplies the flavor-singlet vector current, under which the ∆ ++ baryon is charged,is to cancel this contribution. Therefore, the magnitude of f V is determined by the size ofdisconnected diagrams with the insertion of a single vector current. In perturbation theory, f V is therefore of order g , because at least three gluons must be emitted from the quark loop One-gluon exchange does not contribute because the color factor vanishes. To see that the two-gluon exchangealso vanishes, one may use the γ -hermiticity of the quark propagator, γ S ( x, y ) γ = S ( y, x ) † , the fact that thefree quark propagator S ( x, y ) is Hermitian for fixed ( x, y ) and γ γ µ γ = − γ µ , to show that the two orientationswith which the quark propagators contribute to the quark loop come with opposite signs and cancel each other. (cid:113) y V (13) k O (31)ext = (cid:113) t y A (12)0 A (23) k O (31)ext − (cid:113) t y z A (12)0 A (23) k O (31)ext Figure 1: The chiral Ward identity in the continuum and in the limit m = 0. Continuous horizontallines indicate that the operator is projected onto vanishing spatial momentum. C. Axial Ward identities and determination of c V Once the renormalization factor Z V and improvement coefficients b V and b effV are known,the improvement coefficient c V can be determined by enforcing an axial Ward identity. In thecontinuum theory, the latter can be derived from the specific transformation δψ ( x ) = − α ( x ) γ ψ ( x ) , δψ ( x ) = − ψ ( x ) α ( x ) γ . (28)As the operator to be chirally rotated, we choose O ( y ) = A (23) µ ( y ), and we have δA (23) µ ( y ) = α ( y ) V (13) µ ( y ) . (29)Choosing α ( x ) to be unity inside the slab t < x < t and zero outside, the variation of theaction (per unit α ) is given by δS (12) = − (cid:90) t t d x (cid:90) d x (cid:16) m R, P (12) R ( x ) − ∂ µ A (12) R,µ ( x ) (cid:17) , (30)with t < y < t . We perform the integral over the divergence of the axial current explicitlyin the continuum using Gauss’s theorem. With the additional constraint z / ∈ [ t , t ], we thenobtain the following Ward identity (cid:90) d y (cid:68) δS (12) A (23) R,k ( y , y ) O (31)ext ( z , ) (cid:69) = (cid:90) d y (cid:68) V (13) R,k ( y , y ) O (31)ext ( z , ) (cid:69) , (31)valid in the continuum [4], and impose it to hold on the lattice, at fixed quark mass and barelattice coupling, up to higher order corrections O( a ). The Ward identity in the chiral limitis illustrated in Fig. 1. For each discretization of the vector current ( α = l, cs ), we define ourestimatorˆ c α V ( m q , g ) = (cid:104) (cid:80) y δS (12) A (23) R,k ( y , y ) O (31)ext ( z , ) − ˆ Z (13) V (cid:80) y V α, (13) k ( y , y ) O (31)ext ( z , ) (cid:105) ˆ Z (13) V (cid:104) (cid:80) y a ∂ ν Σ (13) kν ( y , y ) O (31)ext ( z , ) (cid:105) (32)with ˆ Z (13) V ( g , m avq , m q , ) = Z V ( g ) (cid:16) b effV ( g ) am avq + b V ( g ) am q , (cid:17) (33)for the local vector current and ˆ Z (13) V ( g , m avq , m q , ) = 1 for the conserved vector current. InEq. (32) we use the symmetric lattice derivative, as in Eq. (14). Other choices would differ byan O( a ) ambiguity in the definition of c V . The renormalized axial and pseudoscalar operatorsfor i (cid:54) = j are, respectively, given by P ( ij ) R ( x ) = Z P (1 + 3 b P am avq + b P am q ,ij ) P ( ij ) I ( x ) , (34a) A ( ij ) k,R ( x ) = Z A (1 + 3 b A am avq + b A am q ,ij ) A ( ij ) k,I ( x ) , (34b)in terms of the improved operators defined in Eq. (6). The renormalized quark mass is definedthrough the relation m ij = m R,ij Z P (1 + 3 b P am avq + b P am q ,ij ) Z A (1 + 3 b A am avq + b A am q ,ij ) ( i (cid:54) = j ) , (35)such that the combination m R, P (12) R is insensitive to Z P , b P and b P . The renormalizationfactor Z A and the improvement parameters b A and b A have been determined non-perturbativelyin Refs. [10–12].In Eq. (32), the operator O ext can be either the vector O ext = V (31) k ( z , ) or the tensoroperator O ext = Σ (31) k ( z , ) 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 theWard identity vanishes in the chiral limit; on the lattice, the improvement term then involves thetwo-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 seein the next section, in the non-perturbative QCD vacuum, both choices are equally well suitedfor separations between the operators of order of 0.5 fm, because the vector-tensor correlation isthen nonvanishing even in the continuum.There is one subtlety here. In Eq. (30), we sum over all time slices in the range [ t , t ] whichimplies the presence of a contact term for x = y . Therefore, on-shell O( a )-improvement is notsufficient to remove all O( a ) contributions and the limit m → c V for different light quark masses usingEq. (32) and then extrapolating to the chiral limit.Finally, in Appendix A we briefly describe a way to determine the improvement coefficient c A using an axial Ward identity. Our non-perturbative determination of c A , which we can compareto the literature [24], serves as a cross-check of our numerical setup. D. Known perturbative results
The known perturbative results in QCD with N c colors and N f flavors of quarks are thefollowing. The result b g = 0 . N f g + O( g ), independently of the pure gauge action,was obtained in [13]. For degenerate quarks, only the combination b V + N f b V appears, and theperturbative series for b V starts at order g . For the tree-level improved L¨uscher-Weisz action,the results are ( C F = N − N c ) [14, 15] Z V = 1 − . × C F g + O( g ) , (36a) b V = 0 . × C F g + O( g ) , (36b) c l V = − . × C F g + O( g ) . (36c)The tree-level value of c cs V is . III. NUMERICAL SETUP
We use the N f = 2 + 1 CLS (Coordinated Lattice Simulations) lattice ensembles [16] whosemain parameters are given in Table. I. They have been produced using the openQCD code of Ref. [17] using the Wilson-Clover action for fermions and the tree-level Symanzik improvedgauge action. The parameter c SW has been determined non-perturbatively in Ref. [18]. Weconsider four values of the bare coupling β = 3 . , . , .
55 and 3 .
70 which correspond tolattice spacings in the range 0 . − .
085 fm [19]. Ensembles using (anti) periodic boundaryconditions (PBC) and open-boundary conditions (OBC) in the time direction have beengenerated on three different chiral trajectories. Two trajectories with constant m avq and http://luscher.web.cern.ch/luscher/openQCD/ T a b l e I : P a r a m e t e r s o f t h e s i m u l a t i o n s : t h e b a r ec o up li n g β = / g , t h e l a tt i ce r e s o l u t i o n , t h e h o pp i n g p a r a m e t e r κ , t h e p i o n m a ss m π a nd t h e P C A C m a ss . ˆ Z ( ) V i s d e fin e d t h r o u g h E q . ( ) : t h e t w o c o l u m n s c o rr e s p o nd t o t h e li g h t a nd s t r a n g e s p ec t a t o r q u a r k s r e s p ec t i v e l y . T h e l a tt i ce s p a c i n g a i nph y s i c a l un i t s w a s d e t e r m i n e d i n [ ]: a β = . = . f m , a β = . = . f m , a β = . = . f m a nd a β = . = . f m . T r a j ec t o r y C L S β L × T κ l κ s m π [ M e V ] a m P C A C ˆ Z ( ) V ˆ Z ( ) V ˆ c l V ˆ c c s V ˆ c A m a v q = c o n s t . w i t h O B C H . × . . . ( ) . ( ) . ( ) + . ( ) + . ( ) − . ( ) H × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) H × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) N × . . . ( ) . ( ) . ( ) ××× C × . . . ( ) . ( ) . ( ) ××× S . × . . . ( ) . ( ) . ( ) + . ( ) + . ( ) − . ( ) N × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) H . × . . . ( ) . ( ) . ( ) + . ( ) + . ( ) − . ( ) N × . . . ( ) . ( ) . ( ) ××× N × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) N × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) D × . . . ( ) . ( ) . ( ) ××× N . × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) N × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) J × . . . ( ) . ( ) . ( ) − . ( ) + . ( ) − . ( ) m q , s = m ph y s q , s w i t h O B C H . × . . . ( ) . ( ) . ( ) ××× H × . . . ( ) . ( ) . ( ) ××× C × . . . ( ) . ( ) . ( ) ××× N . × . . . ( ) . ( ) . ( ) ××× N × . . . ( ) . ( ) . ( ) ××× D × . . . ( ) . ( ) . ( ) ××× N . × . . . ( ) . ( ) . ( ) ××× m q , l = m q , s w i t h P B C r q c d . × . . . ( ) . ( ) . ( ) ××× B × . . . ( ) . ( ) . ( ) + . ( ) + . ( ) − . ( ) r q c d × . . . ( ) . ( ) . ( ) ××× X × . . . ( ) . ( ) . ( ) ××× X . × . . . ( ) . ( ) . ( ) ××× X × . . . ( ) . ( ) . ( ) ××× m q , l = m q , s w i t h O B C H . × . . . ( ) . ( ) . ( ) ××× H × . . . ( ) . ( ) . ( ) ××× N . × . . . ( ) . ( ) . ( ) ××× N × . . . ( ) . ( ) . ( ) ××× m q ,s = m physq ,s can be used to extrapolate results to the physical limit with physical masses ofthe light and strange quarks. A third trajectory uses degenerate light and strange quarks with m q ,l = m q ,s . Concerning c V , it is enough to consider ensembles on a single chiral trajectory(e.g. m avq = const). However, to determine the two improvement coefficients b V and b effV , wehave to consider at least two different chiral trajectories.For the calculation of the renormalization factor Z V , we need to compute the following three-point correlation function, projected on vanishing momentum C P V P ( x , y ; z ) = a (cid:88) x , y (cid:104) P ( ij ) ( x , x ) V ( jj )0 ( y , y ) P ( ij ) † ( z , ) (cid:105) , (37)and two-point correlation functions C P P ( x , z ) = a (cid:88) x (cid:104) P ( ij ) ( x , x ) P ( ij ) † ( z , ) (cid:105) , (38a) C AP ( x , z ) = a (cid:88) x (cid:104) A ( ij )0 ( x , x ) P ( ij ) † ( z , ) (cid:105) . (38b)Correlation functions are calculated using U(1) stochastic sources with time dilution [20]. Oneach gauge configuration, we generate an ensemble of N s stochastic sources with support on asingle time slice as well as satisfyinglim N s →∞ N s N s (cid:88) s =1 η aα ( x ) s (cid:104) η bβ ( y ) s (cid:105) ∗ = a − δ αβ δ ab δ x,y , (39)where each component is normalized to one, η aα ( x ) ∗ [ r ] η aα ( x ) [ r ] = 1 (no summation). This canbe 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. Inthis case, the source is placed at z = T / t = 0) and the two-pointcorrelation functions are obtained for all values of x ∈ [0 , T / t = T ). For the three-point correlation function, the sink time is placedat x = 3 T / y .For each stochastic source s with support in time slice z , we solve the Dirac equation anddenote the solution vector Φ si ( x ; z ) = a (cid:80) z S ( x, z ) η si ( z ). Correlation functions are given by C ( ij ) P P ( x , z ) = a V (cid:88) x , z Tr (cid:104) S i ( x, z ) † S j ( x, z ) (cid:105) = a N s V (cid:88) s, x Φ si ( x ; z ) † Φ sj ( x ; z ) , (40a) C ( ij ) AP ( x , z ) = a V (cid:88) x , z Tr (cid:104) S i ( x, z ) † γ S j ( x, z ) (cid:105) = a N s V (cid:88) x Φ si ( x ; z ) † γ Φ sj ( x ; z ) , (40b) C ( ij ) P V P ( x , y ; z ) = a V (cid:88) x , y , z Tr (cid:104) S i ( x, z ) † S j ( x, y ) γ S j ( y, z ) (cid:105) (40c)= a N s V (cid:88) x (cid:101) Φ s † ji ( y ; x , z ) γ Φ sj ( x ; z ) , Table II: Values of z , t , t and y for the calculation of the three-point correlation function as definedin Eq. (31). In the last column, we give the two values of y used to interpolate to a line of constantphysics as explained in the text. Note that ensembles at β = 3 .
46 were generated using periodic boundarycondition in the time direction whereas other ensembles were generated using open boundary conditions. β T / (2 a ) z /a [ t , t ] /a y /a − − − − − − . . . . .
775 0 10 20 30 40 50 60 / R ( t , t ) t /a N300 − . − . . . . ˆ c V y /a N300
Figure 2: Plateaus for the ratio R ( t = T , t ) and ˆ c V defined through Eqs. (41) and (32) for the latticeensemble N300. . . . . . . − . − .
002 0 0 .
002 0 .
004 0 . Z V = 0 . b V = 1 . b eff V = 0 . c Z ( ) V am q,l β = 3 . Tr[ M q ] = 0 . M q ] = 0 . M q ] = 0 . M q ] = 0 . . . . . . . − . − .
001 0 0 .
001 0 .
002 0 .
003 0 .
004 0 .
005 0 . Z V = 0 . b V = 1 . b eff V = 0 . c Z ( ) V am q,l β = 3 . Tr[ M q ] = 0 . M q ] = 0 . M q ] = 0 . M q ] = 0 . M q ] = 0 . M q ] = 0 . Figure 3: Results of the fits used to determine the renormalization constant Z V and improvement coeffi-cients b V and b effV using the fit ansatz (42) for two different values of the bare coupling. with V = L the spatial volume. We have used the γ -Hermiticity of the fermionpropagator S ( x, y ) = γ S ( y, x ) † γ and (cid:101) Φ sji is a sequential propagator given by (cid:101) Φ sji ( y ; x , z ) = a (cid:80) x , z γ S j ( y, x ) γ S i ( x, z ) η s ( z ). In practice, since the stochastic sourcesdo 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 pointsource method with a single source location.To compute the correlation functions in Eqs. (30) and (31) we instead use point sourcesand the method of sequential propagators for the three-point correlation functions. A pointsource is first created on time slice z . Then, a sequential inversion is performed using thevariation of the action between time slices t and t as a sequential source. We thereby haveaccess to all y values in the range [ t , t ]. To increase statistics, we also average over equivalentpolarizations k = 1 , ,
3. The values of t , t , z and y used in our simulations are summarizedin Table II. We have computed the correlation functions entering Eq. (31) to leading order inlattice perturbation theory (see Appendix B) in order to test our lattice QCD code. IV. RESULTSA. Results for Z V , b V and b effV Away from the boundary, it is convenient to use the variables t = x − z and t = y − z .For each ensemble, the value of ˆ Z V is estimated from the ratio of three- to two-point correlationfunctions, defined through Eq. (24) with a local vector current. We choose j = l (spectator1quark), which define our renormalization scheme. The ratio has the asymptotic behavior R ( t, t ) −−−−−−−→ t , t − t →∞ Z (12) V , (41)and is fitted to a constant in the plateau region where discretization effects are small. For ensem-bles with antiperiodic boundary conditions in time, we use C P V P ( x , y ; z ) → C P V P ( x , y ; z ) − C P V P ( x , y + T ; z ) to impose the vector Ward identity on each gauge configuration which canhave a nonzero charge due to thermal fluctuations. Typical plots for the ensembles N200 andN300 are given in Fig. 2, and the results for all ensembles are summarized in Table I. In a secondstep, the renormalization constant Z V , and the improvement coefficients b V and b effV at a givenvalue of the bare coupling g are obtained using the fit ansatzˆ Z (12)V ( g , m avq , m q , ) = Z V ( g ) (cid:16) b effV ( g ) am avq + b V ( g ) am q , (cid:17) . (42)The ensembles included in the fit satisfy | am q,l | < .
01 and am avq < .
01 such that higher ordercorrections are expected to be small. The results for each value of β are given in Table III andthe statistical error includes the error on κ cr . The fits for two values of β are shown in Fig. 3.We note that the coefficient b V is significantly larger than the one-loop perturbative estimategiven in Eq. (36b) and that b effV (cid:28) b V . This was expected since the perturbative value startsonly at two loops in perturbation theory. We provide the covariance matrices for the differentvalues of the coupling considered in this workcov( Z V , b V , b effV ; β = 3 .
40) = +6 . × − − . × − − . × − − . × − +1 . × − +1 . × − − . × − +1 . × − +5 . × − , (43a)cov( Z V , b V , b effV ; β = 3 .
46) = +4 . × − +9 . × − − . × − +9 . × − +1 . × − − . × − − . × − − . × − +2 . × − , (43b)cov( Z V , b V , b effV ; β = 3 .
55) = +3 . × − − . × − − . × − − . × − +9 . × − +1 . × − − . × − +1 . × − +1 . × − , (43c)cov( Z V , b V , b effV ; β = 3 .
70) = +3 . × − +2 . × − − . × − +2 . × − +1 . × − − . × − − . × − − . × − +3 . × − . (43d)Finally, we perform a Pad´e fit to obtain the renormalization factor and the improvement coeffi-cients as a function of the bare coupling g . Our final results read Z V ( g ) = 1 − . g × p g + p g p g , (44a) b V ( g ) = 1 + 0 . g × p g p g , (44b) b effV ( g ) = p g p g , (44c)2 . . . .
020 0 .
001 0 .
002 0 .
003 0 .
004 0 .
005 0 .
006 0 .
007 0 . a [fm ]Scaling for Z V Z l V , sub /Z V y ( x ) = 1 − . e a + 1 . e a y ( x ) = 1 − . e a + 5 . e a . . . . .
10 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 .
08 0 . a [fm]Scaling for b V b l V /b FV b l V /b s V . . . .
540 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 .
07 0 .
08 0 . y ( x ) = 1 − . e a + 143 . e a y ( x ) = 1 − . e a + 46 . e a a [fm]Scaling for e b V e b l V / e b FV e b l V / e b s V Figure 4: Left panel: Continuum limit of the ratio Z l V , sub /Z V where Z l V , sub was computed in [12] andwhere Z V refers to our own determination. Middle panel: b l V and b s V correspond to our determinationswith the light or strange spectator quarks, respectively, and b F V is the value computed in [21]. Rightpanel: b F V is the value computed in [21]. We have defined the dimensionless parameter (cid:101) a = a/ . which automatically reproduce the one-loop calculations and where the parameters and covari-ance matrices are given by p ( Z V ) = − . − . − . , cov( Z V ) = +1 . . . . . . . . . × − , (45a) p ( b V ) = (cid:32) − . − . (cid:33) , cov( b V ) = (cid:32) +36 . . . . (cid:33) × − , (45b) p ( b effV ) = (cid:32) +0 . − . (cid:33) , cov( b effV ) = (cid:32) +1 . . . . (cid:33) × − . (45c)To ensure that O( a ) ambiguities vanish smoothly toward the continuum limit, the renormal-ization of the vector current must be performed along a line of constant physics (LCP). Sincethe CLS ensembles have different volumes, we checked explicitly that the impact on the renor-malization factor is extremely small. The observable ˆ Z (12)V has been computed on two sets ofensembles (H105/N101 and H200/N202) generated using the same lattice parameters but withdifferent volumes and the results quoted in Table I are in perfect agreement within statisticalerrors. Second, the correlation functions in Eq. (24) are computed with a source located at z = T / z = T / T / − a , T / − a , and the results areˆ Z (12)V = 0 . . . j = s for the spectator quark, and the values of theimprovement coefficients b V and b effV would differ by an O( a )-ambiguity. To study this effect, wehave repeated the analysis with j = s , and the results are given in Table III. We do not observeany difference for the renormalization constant Z V at our level of precision (both results shoulddiffer only by an O( a ) ambiguity). For b V and b effV , we observe variations by factors of at most1 .
07 and 1 .
7, respectively. In Fig. 4, the continuum limit behavior of the ratio between the tworesults obtained using either a light or a strange spectator quark is shown in blue. For b V , weobserve the expected linear scaling with the lattice spacing. For b V , the ratio goes to one onlyif one includes higher order discretization effects, which appear to be sizable.
1. Comparison of results with previous work
The renormalization factor Z V has been determined independently in [12] using the chirallyrotated Schr¨odinger functional framework. In Fig. 4, we plot the ratio Z l V , sub /Z V where Z l V , sub is extracted from [12] using the line of constant physics called L and where the denominator3 − . − . . . . . .
002 0 .
004 0 .
006 0 .
008 0 .
01 0 . am β = 3 . b c l V b c c V − . − . . . . . .
001 0 .
002 0 .
003 0 .
004 0 .
005 0 .
006 0 . am β = 3 . b c l V b c c V Figure 5: Chiral extrapolations of the improvement coefficients c l V and c c V , respectively for the local andconserved vector currents, for two different values of the bare coupling. corresponds to our own determination. This ratio goes rapidly to one in the continuum limit,even though the expected O( a ) scaling is not observed. However, the maximum deviation,obtained at β = 3 .
40, is less than 1.6%. Empirically, the available data for the departure ofthe ratio from unity can be described by the sum of a linear term and a quadratic term in thelattice spacing (not expected theoretically), or by the sum of a quadratic term and a cubic term.The latter fit in fact describes the data slightly better, see Fig. 4. It also yields coefficients ofreasonable size if one evaluates the lattice spacing say in units of 0.5 fm.The coefficients b V and b effV have also been determined recently in Ref. [21] using a differentsetup, based on the QCD Schr¨odinger functional. A comparison with our results is given inFig 6. For b V , we observe a deviation of about 5%, similar to the O( a ) dependence on thespectator-quark estimated above. In Fig. 4, we show the continuum limit behavior of the ratiowith our own results and we observe the expected linear scaling. However, for b effV , the differencewith the results quoted in Ref. [21] is significant, especially at large couplings g . Again, as canbe 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 froma large ambiguity.From a practical point of view, one should remember that a typical value of am avq is 0.005on the β = 3 .
55 ensembles, so that with 3 b effV (cid:39) .
16 even a 100% ambiguity on b V has animpact below the permille level. Conversely, it could be that O( a ) effects compete with theseterms, resulting in a substantial O( a ) contamination in our determination of b effV . For thephysics applications discussed in the Introduction, it is interesting to compare our values forthe renormalization factor ˆ Z (12)V of the isovector current to those one would obtain using the Z V values of [12] and the b V and b effV values from [21]. We find that our direct estimates are alwaysslightly lower, and that the relative difference depends almost only on g : it is about 1.3% at β = 3 .
40, 0.8% at β = 3 .
46, 0.37% at β = 3 .
55 and 0.12% at β = 3 .
70. We conclude that thesedifferences are of reasonable size, compatible with the expected a (and higher order) ambiguityintroduced by the choice of a specific renormalization condition. B. Results for the improvement coefficient c V For y not too close from t and t , we can extract the value of ˆ c V for each lattice ensemble.In practice, since we want to use a line of constant physics, we choose y − z = 0 .
77 fm andinterpolate linearly between two neighboring time slices when necessary. Deviations from LCPdue to the different sizes of the lattices used should be very mild since we are working in largevolumes. The values of z , t and t in Eqs. (30) and (31) as well as the values of y used in theinterpolation are quoted in Table II. Similar results are obtained using either the vector operatoror the tensor operator as a probe operator in Eq. (31). In practice, we use the linear combination O ext = V (31) k ( z , ) + Σ (31) k ( z , ) which helps to improve the statistical precision. For Z A , we4 . . . . . . . . . . . . . g Z V ( g ) Z V . . . . . . . . . . . g b V ( g ) b V . . . . . . . . . . . . . g b eff V ( g ) b eff V Fritzsch1-loop PT − . . . . . . . . . . . . g c V ( g ) c l V c c V Figure 6: The dependence of the renormalization constant Z V and improvement coefficients b V , b effV and c V on the bare coupling g . The blue and red points correspond to the local and conserved vector currentsrespectively. The plain lines and error bands correspond to our Pad´e fits. For b V and b effV we also compareour results with previous lattice determinations [10, 21]. use the results called Z l A , sub using the L -LCP from [12]. For b A we used the published valuesin [10], and for b A we use values from [22] (in practice, we use b effA which includes the b g -term for Z A , as in Eq. (17)). The results for ˆ c V for each ensemble are given in Table I and differ from c V by the presence of a contact term which vanishes in the chiral limit, as explained in Sec. II C.The improvement coefficient c V is obtained using a linear extrapolation in the light-quark PCACmass m at constant m avq , and the results for each value of the bare coupling are summarizedin Table III. As can be seen in Fig. 5, we observe a very mild chiral dependence. The erroris dominated by the statistical precision of the correlation functions and the error on b A . Theuncertainty on Z A and b A appears to have a negligible impact at our level of precision. Finally,we perform linear or quadratic fits in g to determine c V as a function of the bare coupling. Theresults for the local and the conserved vector currents read c l V ( g ) = − . C F g × (1 + 0 . g ) , (46a) c cs V ( g ) = 12 × (1 − . g ) . (46b)The parametrization is consistent with the perturbative predictions collected in Sec. II D. Ourvalues for c l V are significantly smaller (in magnitude) than the preliminary values determinedin [23] by applying a similar improvement condition in the Schr¨odinger functional. It couldbe due to a large O( a ) ambiguity in the definition of c V . However, our values for Z V alsodiffer by more that one standard deviation from the ones computed in [23]. Since Z V entersin the determination of c V , this could partly explain the disagreement. For example, if we usethe preliminary results for Z V given in [23], our value would be ˆ c l V = − .
146 for H105 andˆ c l V = − .
178 for N200. These observations highlight the need for good control over Z V todetermine precisely c l V . Since the point-split vector current is conserved, Z c V = 1, this issue isabsent for the determination of c cs V .5 V. CONCLUSION
We have determined non-perturbatively the renormalization constant and improvement co-efficients of the local and point-split nonsinglet vector currents with N f = 2 + 1 O( a )-improvedWilson quark action and the tree-level Symanzik improved gauge action. Only one coefficient, f V , is missing but is also expected to be small, as it starts at O( g ) in perturbation theory; inthis regard, we note that for the two other mass-dependent improvement coefficients, b V and b V , the hierarchy expected from perturbation theory is indeed observed in our non-perturbativeresults. All these parameters are required for the full O( a )-improvement of the vector currentand the reduction of discretization effects in lattice simulations. They are essential in the cal-culation of the hadronic vacuum polarization (HVP) contribution to the muon g −
2, where aprecision below 1% is aimed at in the near future.Full O( a )-improvement of the vector current requires one to consider the renormalizationfactor Z V ( (cid:101) g ) at the value of the renormalized coupling (cid:101) g instead of the bare coupling g .We have taken this difference into account by replacing the improvement coefficient b V by theeffective parameter b effV , thus avoiding the use of the unknown coefficient b g .We have obtained the renormalization factor and improvement coefficients by imposing vectorand axial Ward identities at finite lattice spacing and bare quark masses on a set of large volumeensembles. Deviations from the line of constant physics in our renormalization scheme have beenstudied and shown to be small for Z V , b V and b V .Our final results for the different β values used in CLS simulations are summarized in Ta-ble III. We also provide interpolating formulas through Eqs. (44) and (46). As a cross-checkof our methods, we have recomputed the improvement coefficient c A and find good agreementwith the results of Ref. [24], which employ an improvement condition set up in the Schr¨odingerfunctional.Our calculation is the first non-perturbative determination of the improvement coefficients c α V with N f = 2 + 1 Wilson quarks for both the local and (the symmetrized version of) thepoint-split vector currents. The value for the local vector current is small and both values, forthe local and point-split vector currents, are close to their perturbative values.The comparison with the recent findings of Ref. [21] shows that a potentially large O( a )-ambiguity in b V remains, but that it should vanish smoothly in the approach to the continuumlimit. For the vector current renormalization, we find important corrections to the expectedasymptotic O( a ) scaling for the difference between our results and the recent determinationof Ref. [12]. However, we note the relative discrepancy is rather small, and smaller than thatobserved for two different normalization conditions for the axial current [11, 12].In the future, other improvement coefficients may be determined for the lattice action usedhere, thanks to the availability of an extensive set of CLS lattice ensembles. In particular, the N f = 2 + 1 hadronic contribution to the running of the weak mixing angle involves the flavor-singlet vector current, whose improvement coefficient ¯ c V is unknown. A method to determinethe latter based on a chiral Ward identity was proposed in [3]. Acknowledgments
We thank Stefan Sint, Rainer Sommer and Mattia Dalla Brida for helpful discussions andGunnar Bali and Piotr Korcyl for sharing with us preliminary results of the improvement co-efficients of the axial current. We are grateful for the access to the ensembles used here, madeavailable 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 Euro-pean Union’s Horizon 2020 research and innovation programme through Grant Agreement no.771971-SIMDAMA. The CLS lattice ensembles used here were partly generated through com-puting time provided by the Gauss Centre for Supercomputing (GCS) through the John vonNeumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEENat J¨ulich Supercomputing Centre (JSC).6
Table III: Results for the renormalization constant Z V and improvement coefficients b V , b effV and c V fordifferent values of the bare coupling. For Z V , b V and b effV the first (second) line corresponds to the resultsobtained with the light (strange) quark as a spectator quark. For c V , both results for the local andconserved vector currents are given. The value of critical hopping parameter κ cr at β = 3 .
40 is extractedfrom [25]. β κ cr Z V b V b effV b V c l V − . − . − . − . c cs V Appendix A: Determination of the axial improvement coefficient c A In this appendix, we use a similar setup to determine the improvement coefficient c A of theaxial vector current (see Eq. (6)). The latter was previously determined non-perturbativelyin [11] in the framework of the Schr¨odinger functional. This study can be seen as a consistencycheck of our method.Within our numerical setup, described in Sec III, we can replace the axial vector current andthe external operator O ext in Eq. (31) by any other operator without new inversion of the Diracoperator. Therefore, we consider the following axial Ward identity (cid:90) d y (cid:104) δS (12) V (23) R, ( y , y ) O (31)ext ( z , ) (cid:105) = (cid:90) d y (cid:104) A (13) R, ( y , y ) O (31)ext ( z , ) (cid:105) , (A1)with O (31)ext = P (31) . The variation of the action is given by Eq. (30), and, similarly to Eq. (31),with the constraint y / ∈ [ t , t ]. Here, V (23) R, and A (13) R, are the renormalized and O( a )-improvedvector and axial vector currents defined in Eq. (6). If one knows Z V , b V , b effV , b A and b A , then c A can be determined by imposing this equation to hold on the lattice up to O( a ) discretizationeffects, as done in Eq. (32) for c V .Using the same procedure as for c V and with the local vector current, we obtain the resultssummarized in Table I for a subset of the ensembles. Similar to c V , the chiral dependence is verymild and the contribution from the contact term in the left-hand side of Eq. (A1) is removedby taking the limit m →
0. We obtain the results quoted in Table IV. As can be seen onFig. 7, our results are close to the values quoted in [11] obtained using a different method. Wepoint out that, in Eq. (A1), the variation of the action δS (12) was computed with the value of c A published in Ref. [11] such that the two determinations are not strictly independent. Sincethis improvement coefficient has an O( a )-ambiguity, we can attribute this small difference to thedifferent schemes used. Table IV: Results for the improvement coefficient c A at different values of the bare coupling. β c A − . − . − . − . − . − . − . − . − . − . − . . .
02 0 0 .
001 0 .
002 0 .
003 0 .
004 0 .
005 0 .
006 0 . am β = 3 .
55 ˆ c A − . − . − . − . . . . . . . . . g c A ( g ) c A ALPHA [PRD 16]1-loop PT
Figure 7: Left: Extrapolation of c A for one value of the bare coupling g . Right: Improvement coefficient c A as a function of the bare coupling g with a (1,1)-Pad´e model. We also plot the results obtained bythe ALPHA Collaboration in Ref. [11] using a different method. Appendix B: Axial Ward identities in the free theory
In the following, we give the tree-level expressions in lattice perturbation theory for thecorrelation functions involved in the chiral Ward identities, Eqs. (31) and (A1). We have usedthese 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, onthe other hand we restrict ourselves to the equal-mass case, m = m = m = m . At order g with a Wilson quark action, the correlation functions do not depend on c sw . We use thestandard notation ˆ p µ = 2 a sin ap µ , p ◦ µ = 1 a sin ap µ , (B1)as well as A ( p ) = 1 + am + 12 a ˆ p , (B2) B ( p ) = m + (1 + am )ˆ p + 12 a (cid:88) k 2) = B ( p ) A ( p ) , (B5)1 a sinh ( aω p ) = p ◦ + C ( p ) . (B6)The free fermion propagator in the time-momentum representation reads, with the conventionsign(0) = 0, S ( x, y ) ≡ (cid:104) ψ ( x ) ¯ ψ ( y ) (cid:105) = (cid:90) B d p (2 π ) e − ω p | x − y | + i p · ( x − y ) D p (B7) (cid:18) sgn( x − y ) 1 a sinh( aω p ) γ − i γ · p ◦ + C ( p ) + δ x ,y a sinh( aω p ) (cid:19) , where (cid:82) B denotes integration over the Brillouin zone, − πa < p i < πa .8In all three-point functions in this appendix, we assume z < min( y , x ). The correlationfunctions relevant to the Ward identity (A1) are a (cid:88) y e − i p · ( y − z ) (cid:68) A ( y ) P ( z ) (cid:69) (B8)= 4 sign( y − z ) (cid:90) B d q (2 π ) C ( q ) sinh ω p + q + C ( p + q ) sinh( ω q ) D q D p + q e −| y − z | ( ω p + q + ω q ) ,a (cid:88) x , z e i p · ( z − y ) (cid:68) A ( x ) V ( y ) P ( z ) (cid:69) (B9)= 8 (cid:90) B d q (2 π ) e − ω q θ ( x − y )( x − y )+( z − y )( ω p + q + ω q ) D q D p + q f ( p , q ) , with f ( p , q ) = C ( q ) q ◦ · k ◦ − C ( p + q ) q ◦ (cid:12)(cid:12)(cid:12) k = p + q x < y ,C ( q ) (cid:16) sinh( ω q ) sinh( ω p + q ) + q ◦ · k ◦ + C ( q ) C ( p + q ) (cid:17) k = p + q x > y (B10)In Eq. (B9), we have assumed x (cid:54) = y . Finally, under the same assumptions we obtain G ( x − z , y − z , p ) ≡ a (cid:88) x , z e i p · ( z − y ) (cid:104) P ( x ) V ( y ) P ( z ) (cid:105) (B11)= − θ ( x − y ) (cid:90) B d q (2 π ) e − ω q θ ( x − y )( x − y )+( z − y )( ω p + q + ω q ) D q D p + q g ( p , q ) , with g ( p , q ) = sinh( ω q ) (cid:16) sinh( ω q ) sinh( ω k ) + q ◦ · k ◦ + C ( q ) C ( k ) (cid:17) k = p + q . (B12)The special case x = y must be treated separately, G ( y − z , y − z , p ) = − (cid:90) B d q (2 π ) e − ( y − z )( ω p + q + ω q ) D q D p + q (B13) × sinh( ω q )2 (cid:16) [sinh( ω q ) + C ( q )][sinh( ω k ) + C ( k )] + q ◦ · k ◦ (cid:17) k = p + q . The correlation functions relevant to the Ward identity Eq. (31) are a (cid:88) y e i p · ( z − y ) (cid:68) V l, ( y ) Σ ( z ) (cid:69) = − y − z ) (B14) × (cid:90) B d q (2 π ) e − ( ω p + ω p + q ) | z − y | D q D p + q (cid:16) C ( q ) sinh ω p + q + C ( p + q ) sinh ω q (cid:17) , and, for x (cid:54) = y and k ⊥ = ( k , k ), (cid:96) ⊥ = ( (cid:96) , (cid:96) ), a (cid:88) x , z e i p · ( z − y ) (cid:68) A ( x ) A ( y ) Σ ( z ) (cid:69) (B15)= − (cid:90) B d k (2 π ) e − ω k ( x − y ) θ ( x − y ) − ( ω k + ω p − k )( y − z ) D k D p − k f A ( p , k ) ,f A ( p , k ) = 12 (cid:110) sinh( ω k ) (cid:16) C ( k ) sinh( ω (cid:96) ) θ ( x − y ) − sinh( ω k ) C ( (cid:96) ) θ ( y − x ) (cid:17) − C ( k ) k ◦ (cid:96) ◦ + C ( k ) k ◦ ⊥ · (cid:96) ◦ ⊥ + C ( k ) C ( (cid:96) ) (cid:111) (cid:96) = p − k , (B16) f A ( , k ) = (cid:40) C ( k )( k ◦ + C ( k ) ) x > y − C ( k ) k ◦ ⊥ x < y . (B17)9Further, for x (cid:54) = y , we have G A ( x − z , y − z , p ) = a (cid:88) x , z e i p · ( z − y ) (cid:68) P ( x ) A ( y )Σ ( z ) (cid:69) (B18)= − θ ( x − y ) (cid:90) B d k (2 π ) e − ω k ( x − y ) θ ( x − y )+( ω k + ω p − k )( z − y ) D k D p − k g A ( p , k ) ,g A ( p , k ) = − 12 sinh( ω k ) (cid:16) sinh( ω k ) sinh( ω (cid:96) ) − k ◦ (cid:96) ◦ + k ◦ ⊥ · (cid:96) ◦ ⊥ + C ( k ) C ( (cid:96) ) (cid:17) (cid:96) = p − k , (B19) g A ( , k ) = − sinh( ω k )( k ◦ + C ( k ) ) . (B20)Finally, for x = y , G A ( y − z , y − z , p ) = − (cid:90) B d k (2 π ) e ( ω k + ω p − k )( z − y ) D k D p − k sinh ω k (B21) × (cid:104) k ◦ (cid:96) ◦ − k ◦ ⊥ · (cid:96) ◦ ⊥ − C (cid:96) C k − sinh ω k sinh ω (cid:96) − C (cid:96) sinh ω k − C k sinh ω (cid:96) (cid:105) (cid:96) = p − k . [1] S. Aoki et al. , Eur. Phys. J. C , 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 , 034504 (2006)[hep-lat/0511014].[4] M. Guagnelli and R. Sommer, Nucl. Phys. Proc. Suppl. , 886 (1998) [hep-lat/9709088].[5] T. Bakeyev et al. [QCDSF-UKQCD Collaboration], Phys. Lett. B , 197 (2004) [hep-lat/0305014].[6] T. Harris and H. B. Meyer, Phys. Rev. D , 114503 (2015) [arXiv:1506.05248 [hep-lat]].[7] R. Frezzotti et al. [ALPHA Collaboration], JHEP , 048 (2001) [hep-lat/0104014].[8] M. Luscher, S. Sint, R. Sommer and P. Weisz, Nucl. Phys. B , 365 (1996) [hep-lat/9605038].[9] M. Bochicchio, L. Maiani, G. Martinelli, G. C. Rossi and M. Testa, Nucl. Phys. B , 331 (1985).[10] P. Korcyl and G. S. Bali, Phys. Rev. D , no. 1, 014505 (2017) [arXiv:1607.07090 [hep-lat]].[11] J. Bulava, M. Della Morte, J. Heitger and C. Wittemeier, Phys. Rev. D , no. 11, 114513 (2016)[arXiv:1604.05827 [hep-lat]].[12] M. Dalla Brida, T. Korzec, S. Sint and P. Vilaseca, Eur. Phys. J. C , no. 1, 23 (2019)[arXiv:1808.09236 [hep-lat]].[13] S. Sint and R. Sommer, Nucl. Phys. B (1996) 71 [hep-lat/9508012].[14] S. Aoki, R. Frezzotti and P. Weisz, Nucl. Phys. B (1999) 501 [hep-lat/9808007].[15] Y. Taniguchi and A. Ukawa, Phys. Rev. D (1998) 114503 [hep-lat/9806015].[16] M. Bruno et al. , JHEP , 043 (2015) [arXiv:1411.3982 [hep-lat]].[17] M. L¨uscher and S. Schaefer, Comput. Phys. Commun. , 519 (2013) [arXiv:1206.2809 [hep-lat]].[18] J. Bulava and S. Schaefer, Nucl. Phys. B , 188 (2013) [arXiv:1304.7093 [hep-lat]].[19] M. Bruno, T. Korzec and S. Schaefer, Phys. Rev. D , no. 7, 074504 (2017) [arXiv:1608.08900[hep-lat]].[20] J. Foley, K. Jimmy Juge, A. O’Cais, M. Peardon, S. M. Ryan and J. I. Skullerud, Comput. Phys.Commun. , 145 (2005) [hep-lat/0505023].[21] P. Fritzsch, JHEP , 015 (2018) [arXiv:1805.07401 [hep-lat]].[22] P. Korcyl and G. S. Bali, private communication.[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 , 555 (2015) [arXiv:1502.04999 [hep-lat]].[25] G. S. Bali et al. [RQCD Collaboration], Phys. Rev. D94