Renormalising vector currents in lattice QCD using momentum-subtraction schemes

We examine the renormalisation of flavour-diagonal vector currents in lattice QCD with the aim of understanding and quantifying the systematic errors from nonperturbative artefacts associated with the use of intermediate momentum-subtraction schemes. Our study uses the Highly Improved Staggered Quark (HISQ) action on gluon field configurations that include $n_f=2+1+1$ flavours of sea quarks, but our results have applicability to other quark actions. Renormalisation schemes that make use of the exact lattice vector Ward-Takahashi identity for the conserved current also have renormalisation factors, $Z_V$, for nonconserved vector currents that are free of contamination by nonperturbative condensates. We show this by explicit comparison of two such schemes: that of the vector form factor at zero momentum transfer and the RI-SMOM momentum-subtraction scheme. The two determinations of $Z_V$ differ only by discretisation effects (for any value of momentum-transfer in the RI-SMOM case). The RI$^{\prime}$-MOM scheme, although widely used, does not share this property. We show that $Z_V$ determined in the standard way in this scheme has $\mathcal{O}(1\%)$ nonperturbative contamination that limits its accuracy. Instead we define an RI$^{\prime}$-MOM $Z_V$ from a ratio of local to conserved vector current vertex functions and show that this $Z_V$ is a safe one to use in lattice QCD calculations. We also perform a first study of vector current renormalisation with the inclusion of quenched QED effects on the lattice using the RI-SMOM scheme.


I. INTRODUCTION
Lattice QCD is the method of choice for the accurate calculation of hadronic matrix elements needed for a huge range of precision particle physics phenomenology aimed at uncovering new physics. Compelling evidence of new physics in the comparison of experiment to the Standard Model (SM) has so far proved elusive, however, and this is driving the need for smaller and smaller uncertainties on both sides. This means that the error bars from lattice QCD calculations must be reduced to sub-1% levels.
Here we address uncertainties coming from the renormalisation of lattice QCD operators to match their continuum QCD counterparts. This renormalisation is needed so that the hadronic matrix elements of the operators calculated in lattice QCD can be used in continuum phenomenology. Ideally the uncertainty from the renormalisation factors, Z, should be much less than other lattice QCD uncertainties (such as statistical errors) in the hadronic matrix element calculation.
Defining QCD on a space-time lattice provides an ultraviolet cutoff on the theory of π/a where a is the lattice spacing. This is a different regularisation than that used in continuum formulations of QCD and hence we expect a finite renormalisation to be required to match lattice QCD and continuum operators. This renormalisation takes account of the differing ultraviolet behaviour * d.hatton.1@research.gla.ac.uk † christine.davies@glasgow.ac.uk ‡ http://www.physics.gla.ac.uk/HPQCD in the two cases and hence can be calculated as a perturbative series in the strong coupling constant, α s , at a scale related to the ultraviolet cutoff. Lattice QCD perturbation theory is notoriously difficult, however, and very few renormalisation constants have been calculated beyond O(α s ) (for an example of a two-loop renormalisation in lattice QCD perturbation theory see [1]). It therefore seems clear that this route will not give accurate enough results for the future. Instead we concentrate here on other approaches that can be implemented using results from within the nonperturbative lattice QCD calculation. These approaches will typically still need to make use of perturbation theory to provide a full matching to a preferred continuum scheme such as MS, but if this perturbation theory can be done in the continuum to high order then much improved accuracy should be possible.
At the heart of these nonperturbative-on-the-lattice approaches is always the idea that we can construct a short-distance operator on the lattice whose leading term in an operator product expansion is the operator that we wish to study. The matrix elements that we calculate on the lattice, and use to determine Z, will be dominated by those from the leading operator. There will inevitably be contamination, however, from subleading terms in the expansion, i.e. higher-dimension operators multiplied by inverse powers of some scale. This means then that nonperturbative artefacts can enter the determination of Z and these must be understood and controlled in order to make use of the Z obtained [2].
Here we will study the renormalisation factor Z V associated with the flavour-diagonal vector current that couples to the photon. This current is conserved in contin-uum QCD and has no anomalous dimension. Hence we can study the lattice QCD determination of Z V directly, and its dependence on the lattice spacing, without having to combine it with a matrix element for the vector current determined in lattice QCD. Z V is a special case of a renormalisation constant that can be calculated exactly in lattice QCD, i.e. without the need for any continuum perturbation theory and without nonperturbative artefact contamination. It is important to use a method that allows for such a calculation if we want an accurate normalisation.
It is possible to write down conserved vector currents in lattice QCD and use these, knowing that they do not require renormalisation because there is an exact vector Ward-Takashashi identity. Conserved vector currents are not generally used, however, because they are complicated objects, especially for discretisations of QCD that are highly improved. The removal of tree-level discretisation errors at O(a 2 ) from the covariant derivative in the Dirac equation requires the addition of operators that extend over three links [3]. The conserved current then contains both one-link and three-link terms and this is the case for the Highly Improved Staggered Quark (HISQ) action that we will use here (see Appendix A). We demonstrate explicitly how the vector Ward-Takahashi identity works in this case.
The HISQ action was designed [4] to have very small discretisation effects and this allows its use to study both light and heavy quark phenomenology [5]. Whenever a vector current is needed for phenomenology, however, it is much easier to use a nonconserved local (or simple one-link point-split) vector current than the conserved one [6][7][8]. This must then be renormalised.
Renormalisation schemes for nonconserved currents that make use (not necessarily explicitly) of ratios of matrix elements for conserved and nonconserved vector currents have a special status because nonperturbative contributions from higher dimension operators are suppressed by powers of a 2 . They give renormalisation constants, Z V , for nonconserved lattice vector currents that are exact in the a → 0 limit. Such a Z V can then be combined with a matrix element of that nonconserved current in the lattice QCD calculation and the result extrapolated to zero lattice spacing. The same answer will be obtained in that limit with any such Z V .
Following the discussion of perturbative matching earlier we can think of an exact Z V as consisting of a perturbative series in α s that depends on the form of the vector current (and also in principle terms arising from small instantons [9] or other nonperturbative effects of this kind) plus discretisation effects that depend on the scheme and vanish as a → 0 [10]. Note that we do not need to know what the perturbative series is; the method is completely nonperturbative. Which exact Z V to use is then simply an issue of numerical cost to achieve a given uncertainty and/or convenience.
One standard exact method for renormalising nonconserved vector currents in lattice QCD is to require (elec-tric) charge conservation i.e. that the vector form factor between identical hadrons at zero momentum transfer should have value 1. Since this result would be obtained for the conserved current, Z V is implicitly a ratio of nonconserved to conserved current matrix elements between the two hadrons. This method is numerically fairly costly because it requires the calculation of two-point and threepoint correlation functions. It can give numerically accurate results (O(0.1%) uncertainties) when averaged over a sufficiently large sample (hundreds) of gluon field configurations. As above, we expect the Z V determined from this method (which we will denote Z V (F(0))) to be equal to a perturbative matching factor up to discretisation effects. This was tested by HPQCD in Appendix B of [11] for the local vector current made of HISQ quarks. Values for Z loc V (F(0)) were calculated at multiple values of the lattice spacing and gave a good fit to a perturbative expansion in α s plus discretisation effects, constraining the O(α s ) coefficient to have the known value determined in lattice QCD perturbation theory.
Alternative methods of determining renormalisation factors by defining a variety of momentum-subtraction schemes on the lattice [12][13][14][15] can produce precise results for Z factors at lower computational cost. However, only some of these schemes are exact for Z V in the sense defined above.
The momentum-subtraction schemes define Z V from the ratio of two matrix elements calculated between external quark states of large virtuality, µ 2 , in a fixed gauge. Working at large µ 2 is part of the definition of these schemes because nonperturbative contributions from higher-dimension operators will in general be suppressed by powers of µ 2 and not a 2 as above. A wavefunction renormalisation factor is determined from the quark propagator. A vertex renormalisation factor comes from an amputated vertex function for the vector current, on which momentum-subtraction renormalisation conditions have been imposed. Z V is then obtained as the ratio of these two factors, with tiny statistical errors from a handful of gluon field configurations if 'momentum sources' are used [16].
The momentum-subtraction scheme known as RI-SMOM [15] is constructed around the Ward-Takahashi identity and so designed to give Z V = 1 for the lattice conserved current. We show explicitly that this is true for the HISQ action. This means that implementing the RI-SMOM scheme for nonconserved currents is equivalent to taking a ratio of vector vertex functions for conserved and nonconserved currents. We compare the Z V values obtained in the RI-SMOM scheme, Z V (SMOM), to those from the form factor method for the local vector HISQ current. We are able to show that, as expected, Z loc V (SMOM) differs from Z loc V (F (0)) only by discretisation effects so that the two methods will give the same answer for physical matrix elements in the continuum limit.
A popular momentum-subtraction scheme that does not make use of the vector Ward-Takahashi identity is the RI -MOM scheme [12,13]. We show that in this scheme the Z V values for both the conserved and local vector currents are not exact but have contamination from nonperturbative (condensate) artefacts that survive the continuum limit. To make use of this scheme Z V must be redefined to use instead a ratio of the vector vertex function for conserved and nonconserved currents. We show the results from implementing this method.
We stress here that we are determining Z V very precisely and hence comparing values with uncertainties at the 0.1% level. Previous work has compared values for Z V for nonconserved currents from methods that use Ward identities and the RI -MOM scheme (for example [17,18]) and concluded that there was agreement at the 1% level. Our more accurate results show clear disagreement, most obviously in the analysis for the conserved current.
Our earlier argument that 0.1% accuracy is needed for renormalisation constants in pure lattice QCD can be extended when we study the impact of adding QED effects. When we allow the valence quarks to have electric charge (i.e. adding quenched QED to lattice QCD) we see a tiny impact (less than 0.1%) on Z V using the HISQ action. We can now quantify and analyse this effect using the RI-SMOM scheme, having established that the nonperturbative Z V values behave correctly.
The paper is laid out as follows: We first discuss in Section II the exact lattice vector Ward-Takahashi identity that gives the conserved vector current for the HISQ action. We then give a brief overview of the momentumsubtraction schemes, called RI-SMOM and RI -MOM, that we will use (abbreviating the names to SMOM and MOM) in Section III and, following that, a brief description of our lattice set-up in Section IV. We show how the Ward-Takashashi identity works for the HISQ action in Section IV A so that the conserved current is not renormalised. This is then translated into the RI-SMOM scheme in Section IV B where Z V = 1 is obtained for the conserved current at all a and µ values. For RI -MOM, however, condensate contributions are clearly evident in the Z V values for the conserved current as shown in Section IV C. In Sections IV D and IV E we demonstrate the impact of the protection from the Ward-Takahashi identity on the renormalisation factors for the simple local vector current (ψγ µ ψ with the fields at the same spacetime point). The difference between Z loc V (SMOM) and Z loc V (F(0)) is purely a discretisation effect; in the RI -MOM scheme we demonstrate how to achieve the same outcome with a renormalisation factor that is a ratio between that for the local and conserved currents. In Section V we show the impact of quenched QED on the Z V values obtained for the local current in the RI-SMOM scheme and compare to our expectations based on the work in earlier sections. Finally, in Section VI we discuss the implications of these results for ongoing and future calculations and give our conclusions. A similar picture to that for the local current is seen for the 1-link pointsplit vector current and we give the RI-SMOM results for this case in Appendix B. We reiterate the shorthand notation that we will use for the renormalisation constants for clarity. Z x V (A) renormalises the lattice vector current x (cons, loc, 1link) to match the continuum current (in e.g. MS) and has been calculated in the scheme A (F(0), SMOM, MOM).

II. THE VECTOR WARD-TAKAHASHI IDENTITY ON THE LATTICE
For both continuum and lattice [19,20] actions the derivation of the vector Ward-Takahashi identity proceeds from the observation that the path integral is invariant under a local change of the fermion field variables ψ and ψ (only) that has unit Jacobian. Then An example of such a transformation is to multiply ψ, say at point x, by phase e i and ψ(x) by e −i : Expanding Eq. (1) to first order in and denoting ∆X = X − X to this order gives If we consider the path integral for the two point correlator ψ(y 1 )ψ(y 2 ) then ∆f becomes the difference of propagators from the points y 1 and y 2 to x. ∆S can be recast into the form ∆ µ J µ , allowing us to identify the conserved current J associated with S. We have The right-hand side is zero unless y 1 or y 2 overlaps with x (and not with each other). Note that ∆ µ J µ (x) is centred on the point x.
On the lattice ∆ µ can be a simple forward (∆ + µ ) or backward (∆ − µ ) finite difference over one link. The current J µ must then be chosen appropriately so that We give J + µ for the HISQ action [4] that we use in Appendix A. As discussed in the Introduction, it is rather complicated. It contains a number of 3-link terms because of the Naik term [3] that removes tree-level a 2 errors in the action. The position-space Ward-Takahashi identity of Eq. (4) provides a test of the implementation of the conserved current and we have checked that this works for our implementation exactly on a single gluonfield configuration for a variety of choices of y 1 and y 2 .
We can perform the exact Fourier transform on the lattice of Eq. (4). The left-hand side becomes where a is the lattice spacing and we take q = p 1 − p 2 .
x is the mid-point of the link between x and x +μ. The right-hand side becomes where S is the quark propagator. Then, multiplying on both sides by the product of inverse quark propagators we reach the lattice version of the standard expression for the Ward-Takahashi identity, is the amputated vertex function for the vector current J µ,+ (absorbing a factor of e iaqµ/2 into the vertex function since J µ,+ sits on a link). This equation is exact, gluon-field configuration by configuration, in lattice QCD and we will demonstrate this for the HISQ action in Section IV A.
As is well-known, Eq. (8) tells us that any rescaling of the vertex by renormalisation on the left-hand side has to match rescaling of the inverse propagators on the righthand side. This means that J µ,+ is not renormalised, i.e. that the renormalisation factor for this conserved current, Z cons V =1. Since this is also true for the conserved current in the continuum MS scheme then the matrix elements of the lattice conserved current will agree in the continuum limit with those in the MS scheme.
A renormalised nonconserved vector current, written for example as Z loc V V loc,µ for a local current, obeys the same equations as J µ,+ since it is by definition the same operator up to discretisation effects on the lattice [10]. For the HISQ action Again this is well-known, but we point it out here because it has implications for the accuracy of the determination of Z loc V on the lattice. It means that, if Z loc V is determined by a procedure which uses the Ward-Takahashi identity and gives 1 for the renormalisation of J µ,± , then Z loc V must be free of systematic errors from nonperturbative (condensate) artefacts in the continuum limit because these must cancel between the left-and righthand sides of Eq. 8. Z loc V can in principle be determined by substituting Z loc V V loc into the lefthand side of Eq. (8) for any p 1 and p 2 . Hadronic matrix elements of Z loc V V loc will then differ from the results in the continuum purely by discretisation effects (which will depend on p 1 and p 2 ) that can be extrapolated away straightforwardly using results at multiple values of the lattice spacing. The Z V so obtained is completely nonperturbative.
Using Eq. (8) in its full generality is unnecessarily complicated and there are lattice QCD methods that make use of it in specific, and simpler, kinematic configurations.
As q → 0 the identity of Eq. (8) can be used to show that the vector form factor for the conserved current between quark or hadron states of the same momentum will be unity. The inverse of the vector form factor at the same kinematic point for a nonconserved current then gives its Z V value. This method clearly satisfies the criteria above for an exact determination of Z V .
We now discuss momentum-subtraction renormalisation schemes on the lattice and the extent to which they make use of Eq. (8).

III. MOMENTUM-SUBTRACTION SCHEMES USED ON THE LATTICE
Momentum-subtraction schemes are useful intermediate schemes between the lattice regularisation and the continuum MS scheme in which it is now standard to quote results for scheme-dependent quantities. If the same momentum-subtraction scheme is implemented both in lattice QCD and in continuum QCD then the continuum limit of the lattice results will be in the continuum momentum-subtraction scheme (and should be independent of lattice action at that point). They can then be converted to the MS scheme using continuum QCD perturbation theory.
A momentum-subtraction scheme imposes renormalisation conditions on matrix elements between (in the cases we consider) external quark states so that the treelevel result is obtained, i.e. Z Γ is defined by for some operator O Γ = ψΓψ, and p 1 | and |p 2 external quark states with momenta p 1 and p 2 , typically taken to have large magnitude. This calculation must of course be done in a fixed gauge, and this is usually taken to be Landau gauge, which can be straightforwardly implemented in lattice QCD. Effects from the existence of Gribov copies under the gauge-fixing could arise in general; here we show that there are no such effects for Z V determined using the Ward-Takashahi identity.
Here we will concentrate on the RI-SMOM scheme [14,15]. This scheme uses a symmetric kinematic configuration with only one scale so that p 2 1 = p 2 2 = q 2 = µ 2 (where q = p 1 − p 2 ). The wavefunction renormalisation is defined (using continuum notation) by The vector current renormalisation follows from requiring The traces here are over spin and colour and normalisations are chosen so that Z q = Z V = 1 at tree-level. The equations above are given for the continuum SMOM scheme. On the lattice we must take care to define the appropriate discretisation for q µ and q 2 in the various places that they appear. Below we will see what form q µ must take in combination with the vertex function for the conserved current. The RI-SMOM scheme was defined with the vector Ward-Takahashi identity in mind [15]. This reference shows how the identity defines the projectors needed for the vector vertex function in the continuum (given in in Eq. (12)) so that Z V = 1 for the conserved current. Here we repeat this exercise, but now on the lattice. Returning to the Ward-Takahashi identity in Eq. (8) we can multiply both sides by / q and take the trace whilst dividing bŷ q 2 (withq a discretisation of q to be defined later). This gives We can simplify the right-hand side assuming that the inverse propagator takes the general form S −1 (p) = i / pΣ V (p 2 ) + Σ S (p 2 ) in the continuum (from relativistic invariance). Then, for the SMOM kinematics, On the lattice this formula could be broken by discretisation effects. We do not see noticeable effects of this kind with the HISQ action, however, as we will discuss in Section IV B. Using Eq. (14) in Eq. (13) and multiplying by i then gives, from the Ward-Takahashi identity (15) From Eq. (11) we see that the righthand-side of this expression is Z q in the RI-SMOM scheme. Comparing the left-hand side to Eq. (12) we see that this is Z q /Z cons V in the RI-SMOM scheme where Z cons V is the Z V factor for the conserved current and the Ward-Takahashi identity requires us to discretise q µ as 2 sin(aq µ /2)/a (q is defined in Eq. (21)). Then, from Eq. (15), we expect that Z cons V (SMOM) = 1 on the lattice and no further renormalisation is needed to match to MS. Notice that this works for any value of q.
We will show by explicit calculation that Z cons V (SMOM) = 1 for the HISQ action in Section IV B. This is not true configuration by configuration, however. It does require an average over gluon fields.
Another popular momentum-subtraction scheme is RI -MOM [12,13], abbreviated here to MOM. In this scheme Z q is defined in the same way, by Eq. (11), but Z V is defined by a different projector for the vector vertex function and the kinematic configuration for the MOM case is p 1 = p 2 = p so that q = 0. Instead of Eq. (12) we have, in the MOM scheme, Since this scheme does not correspond to a Ward-Takahashi identity, Z V determined this way needs further renormalisation to match to the MS scheme. More problematically, as we will show in Section IV, Z cons V (MOM) for the HISQ action is significantly different from 1 and is contaminated by nonperturbative condensate effects.
The RI-SMOM γµ scheme [15] is similar to RI -MOM above but uses the SMOM kinematics with p 2 1 = p 2 2 = q 2 . To calculate the renormalisation constants for nonconserved currents we must combine the calculation of the vector vertex function for that current (Eq. (12) and appropriate modifications of it as described in the text) with the calculation of the wave-function renormalisation (Eq. (11)). We describe the results for the HISQ local vector current in the SMOM scheme in Section IV D. We are able to show that the renormalisation factor for the local vector current in the SMOM scheme differs from that using the form factor method purely by discretisation effects, demonstrating that it is an exact form of Z V . The discretisation effects depend on q but the method is exact for any q; this is in contrast to the usual idea of a 'window' of q values to be used in momentum-subtraction schemes on the lattice [12].
The RI -MOM scheme is not exact, as discussed above. We show in Section IV E that a modification of the method (reverting to one of the original suggestions in [12]) does, however, give an exact Z V .
There are technical issues associated with implementing momentum-subtraction schemes for staggered quarks that we will not discuss here. We use the techniques developed in [21] and summarised again in [2] in the context of the RI-SMOM scheme. We will only discuss here specific issues that arise in the context of the vector current renormalisation.

IV. THE LATTICE QCD CALCULATION
We perform calculations on n f = 2 + 1 + 1 gluon field configurations generated by the MILC collaboration [22,23] listed in Table I. These ensembles use an improved gluon action which removes discretisation errors through O(α s a 2 ) [24]. They include the effect of u/d, s and c quarks in the sea using the HISQ action [4].
All gauge field configurations used are numerically fixed to Landau gauge by maximising the trace over the TABLE I. Simulation parameters for the MILC gluon field ensembles that we use, labelled by set number in the first column. β = 10/g 2 is the bare QCD coupling and w0/a gives the lattice spacing [25], using w0 = 0.1715(9) fm [26] determined from the π meson decay constant, fπ. Note that, for each group of ensembles at a given value of β we use the w0/a value corresponding to the physical sea quark mass limit [2], using results from [27]. Ls and Lt give the lattice dimensions. am sea l , am sea s and am sea c give the sea quark masses in lattice units. Set 1 will be referred to in the text as 'very coarse', sets 2-5 as 'coarse', set 6 as 'fine' and set 7 as 'superfine'. This is enough to remove the difficulties related to loose gauge fixing discussed in [2]. We use broadly the same calculational set up as in [2] but here we are considering vector current vertex functions rather than scalar ones. To implement momentumsubtraction schemes for staggered quarks we need to use momenta within a reduced Brillouin zone [21] −π/2 ≤ ap ≤ π/2.
For each momentum ap 1 or ap 2 we then calculate propagators or vertex functions with 16 copies of that momentum, ap 1 + πA and ap 2 + πB where A and B are four-vectors composed of 0s and 1s. This then enables us to do the traces over spin for specific 'tastes' of vector current implied by equations such as Eq. (15). There is also a trace over colour in this equation so the S −1 (q) factor on the righthand side, for example, is actually a 48 × 48 matrix. Where necessary we will use the notation of [21] to denote specific spin-tastes. As an example γ µ ⊗ I is the 16 × 16 matrix of 0s and 1s that projects onto a taste-singlet vector in AB space. Twisted boundary conditions are utilised to give the incoming and outgoing quarks arbitrary momenta [28,29]. For the SMOM kinematics we take, with ordering (x, y, z, t), For the MOM kinematics we take ap 2 = ap 1 . A range of aµ values are chosen at each lattice spacing, satisfying Eq. (17). This allows us to reach µ values of 3 GeV on coarse lattices and 4 GeV on fine and superfine lattices [2]. The µ values can be tuned very accurately (to 3 dec. places).
Relatively small samples (20 configurations) give small statistical uncertainties for Z V at the µ values that we use (with momentum sources for the propagators). A bootstrap method is used to estimate all uncertainties and include correlations between results at different µ values on a given ensemble. Bootstrap samples are formed for each Z q and each Λ V and the bootstrap averages are then fed into the ratio to determine Z V .
All of our results are determined at small but non-zero valence quark mass. Degenerate masses are used for the incoming and outgoing quarks (but note that there is no need for the calculation of disconnected contributions). As the momentum-subtraction schemes that we consider are in principle defined at zero valence quark mass (but direct calculation at this point will have finite-volume issues) it is necessary to calculate each Z V at different quark masses and then extrapolate to the am val = 0 point. To do this we perform all calculations at three masses corresponding to the light sea quark mass on a given ensemble, am l , and at 2am l and 3am l . Dependence on am val can come from discretisation effects and from the contribution of nonperturbative condensate terms.
We follow the procedure used for Z m in [2] and extrapolate Z V results using a polynomial in am val /am s : We find no need for higher powers of am val /am s here as the valence mass dependence of Z V is observed to be very mild in all cases. For the priors for the coefficients d i we use {0 ± 0.1, 0 ± 0.01} at µ = 2 GeV with the widths decreased according to µ −2 .
Any sea quark mass dependence should be suppressed relative to the valence mass dependence by powers of α s and this was observed in [2]. As the valence mass dependence is already negligible the sea mass dependence should be tiny here and we ignore it.

A. The Ward-Takahashi identity on the lattice
In this Section we test the exact lattice Ward-Takahashi identity for HISQ quarks, i.e. Eq. (8). If we have correctly implemented the lattice conserved vector current, this equation is true as a 3 × 3 matrix in colour space. It is also true for any p 1 and p 2 (except that it reduces to 0=0 for p 1 = p 2 ), any values of the quark mass and any gauge. We test it for the SMOM kinematic configuration of Eq. (18). Figure 1 shows the results as a ratio of the difference of inverse propagators on the righthand side of Eq. (8) to the amputated vertex function for the conserved vector current on the lefthand side. This is averaged over colour components (which all agree) and summed over the two non-zero components of q µ (which take the same value  aµ/ √ 2 in each of the y and z directions for the SMOM kinematics). The Ward-Takahashi identity (Eq. (8)) requires this ratio to be exactly equal to 2 sin[aµ/(2 √ 2)], which is plotted as the line.
The plot shows that this expectation works to high precision (double precision accuracy here), on a single configuration taken as an example from Set 2. Results are given for three different aµ values with two different valence quark masses and in two different gauges. The agreement between the points and the line demonstrates the Ward-Takashi identity working explicitly on the lattice for the conserved HISQ current of Eq. (A1). The agreement seen in two different gauges is evidence that the Ward-Takahashi identity works in any gauge, as it must, and therefore its operation is also independent of any Gribov copy issue in the gauge-fixing procedure.

B. ZV for the conserved current in the RI-SMOM scheme
To determine Z V for the HISQ conserved current in the RI-SMOM scheme we adapt Eqs (11) and (12) to the case of staggered quarks on the lattice, as partly discussed already in Section III. For staggered quarks the inverse propagator is a taste-singlet [21] and so the HISQ version  of Eq. (11) is The trace is now over colour and AB-space index.q is given by This choice is dictated by the momentum-subtraction requirement that Z q should be 1 in the non-interacting (tree-level) case and the fact that the derivatives in the HISQ action are improved through O(a 2 ) [4]. Likewise the HISQ calculation for Z V for this case is given by In Sec. III it was shown how the Ward-Takahashi identity leads to the exact expression of Eq. (13) on the lattice when the conserved current is used in the vertex function. In order to obtain Z V = 1 for the conserved current we also need Eq. (14) to be satisfied exactly. In Fig. 2 we give a test of this relationship. The figure shows the ratio of the difference of the two inverse propagators with momentum p 1 and p 2 to that of the propagator with momentum q, where the inverse propagators are multiplied by / q and the trace taken. We useq here (Eq. (21)) instead of simply q to be consistent with what we use in the determination of Z q in Eq. (20) above. The results for the ratio plotted would be the same for q as forq. The results for the ratio in Fig. 2 are seen to be consistent with 1.0 to better than 0.05%. The statistical uncertainties plotted are from a bootstrap over results from 20 gluon field configurations. Figure 2 shows that discretisation effects in the HISQ action have no effect on Eq. (14) at the level of accuracy to which we are working. There are no tree-level a 2 errors with HISQ [4] and there is a U(1) axial symmetry; both of these constrain the form that discretisation effects can take [21]. A further constraint comes from the form for the p 1 and p 2 momenta (and q) to achieve the SMOM kinematics. Each has only two non-zero momentum components, as shown in Eq. (18). This means, for example, that discretisation errors in S −1 containing 3 different γ matrices and associated momenta are zero. Figure 3 shows the resulting Z V value obtained for the conserved vector current in the RI-SMOM scheme, combining the results from Eqs (22) and (20) and performing the extrapolation to zero quark mass as described in Sec. IV (this has very little impact). The value obtained for Z V for the conserved current is 1 to better than 0.05% at all µ values. Fitting the results shown in Fig. 3   The kinematic conditions in the MOM scheme are that the incoming and outgoing quark fields for the vertex function should have the same momentum, ap 1 = ap 2 so that aq = 0. We will denote this momentum by ap with |ap| = aµ. We take the form of ap to be that of ap 1 in the SMOM scheme (Eq. (18)). To implement the RI -MOM scheme we determine the wave-function renormalisation, Z q (p), in the same way as for the RI-SMOM scheme using Eq. (20). To determine Z V we use This uses the RI -MOM vector vertex projector, which is simply γ µ (see Eq. (16)), expressed here in the appropriate taste-singlet form for implementation with staggered quarks. Because the conserved current is a pointsplit operator the tree-level vertex function is not simply 1. We therefore need to divide by the tree-level matrix element for the conserved current that we denote here V cons γµ⊗I . How to calculate these tree-level factors is discussed in [21]. We have The spin-taste 4-vector S − T is composed of 1s and 0s. For the taste-singlet vector it takes value 1 for component µ and 0 otherwise. So the only components of the product that do not take value 1 are those for component µ that matches the direction of the current, provided that ap has a non-zero component in that direction.
Because the RI -MOM scheme is not based on the Ward-Takahashi identity Z V will not be 1 for the conserved current. This means that to reach the MS scheme, even for the continuum RI -MOM scheme, requires an additional renormalisation factor. The renormalisation factor that takes the lattice vector current to the continuum is then Z MOM,raw V is the raw renormalisation factor calculated using Eq. (23) on the lattice. The factor Z MS/MOM V can be determined from the perturbative QCD expansions in the continuum for the conversion between RI -MOM and RI-MOM given in [13] (see [30] and the Appendix of [21]). The values needed for our µ values are given in Table II; they are all close to 1 since the expansion starts at O(α 2 s ). Figure 4 shows our results for Z V for the conserved HISQ current obtained by implementing the RI -MOM scheme on the lattice. We have converted the Z V to the value that takes the lattice results to the MS scheme using Eq. (25). Results are shown, after extrapolation to zero valence quark mass, at a variety of µ values from 2 GeV to 4 GeV and at three different values of the lattice spacing. It is immediately clear that the values of Z cons V (MOM) are not 1. This is in sharp contrast to results in the RI-SMOM scheme where, as we showed in Section IV B, the value 1 is obtained. This result is shown by the black line at 1 in Fig. 4.
To understand the discrepancy from 1 for Z cons V in the RI -MOM case, we fit the points shown in Fig. 4 (including the correlations between them) to a form that allows for both discretisation effects and condensate contributions: cond,a 2 (aΛ/π) 2 ] .
Note that this constrains Z cons V (MOM) to be 1 in the continuum once condensates are removed. Here α MS (µ) is the value of the strong coupling constant in the MS scheme at the scale µ calculated from running the value obtained in [27] using the four-loop QCD β function. The fit allows for discretisation errors of the generic form (aµ) 2i and terms O(α s (aµ) 2i ); only even powers of a appear due to the remnant chiral symmetry of staggered quarks. Note that in principle we have removed (aµ) 2 terms by dividing by V cons γ⊗I ; the fit returns only a small coefficient for this term. The α s -suppressed discretisation terms are included as the very small statistical uncertainties on the results mean that these terms can have an effect in the fit. The fourth term allows for systematic uncertainty from the missing α 4 s term in the RI -MOM to MS conversion factor (Eq. 25).
The condensate terms on the final line of Eq. (26) start at 1/µ 2 to allow for the gauge-noninvariant A 2 condensate present in the operator product expansion (OPE) of the quark propagator [31]. For the MOM kinematic setup it is not possible to perform an OPE for the vertex functions as they are not short-distance quantities (q = 0), so a complete analysis of what nonperturbative artefacts we expect to see in Z V is not possible. However, on general grounds we expect terms with inverse powers of µ to appear and we allow these terms also to have discretisation effects. We include even inverse powers of µ up to 1/µ 8 .
We use a Bayesian fit approach [32] in which coefficients are constrained by priors with a Gaussian distribution of a given central value and width. All coefficients in the fit form of Eq. (26) are given priors of 0 ± 1, except for that of the (α s /π) 4 term which has prior 0 ± 5 based on the lower-order coefficients. The choices for the priors are based on reasonable values for the coefficients of the terms in the fit. For example, discretisation effects are expected to appear as even powers of a physical scale (such as µ or Λ here) divided by the ultraviolet cutoff (π/a) with coefficients of order one.
The results of the fit are shown as the coloured dashed lines in Fig. 4. The fit has a χ 2 /dof of 0.6. It is already obvious from the figure that discretisation effects are not the only source of the discrepancy in Z V from 1. This is emphasised by attempting the fit without condensate terms (i.e. missing the last line of Eq. (26)). Without the condensate terms the quality of the fit is very poor, with a χ 2 /dof of 7.7, in contrast to the fit of Eq. (26). The sizeable contribution from the lowest order condensate is reflected in the coefficient found by the fit of The higher-order condensates cannot be pinned down by the fit. The correct answer for Z V for the conserved current in the continuum limit is, of course, 1. Our results and fit show that this can only be obtained from a calculation in the RI -MOM scheme by working at multiple µ values at multiple values of the lattice spacing and fitting as a function of µ and a to identify and remove the condensate contributions. If this is not done, systematic errors of O(1%) (depending on the µ value) are present in Z V , as is clear from Fig. 4.
The issue will resurface when we discuss the use of the RI -MOM scheme to renormalise nonconserved currents, specifically the HISQ local vector current, in Section IV E.

D. ZV for the local current in the RI-SMOM scheme
We now turn to the calculation of the renormalisation constant for a nonconserved vector current using the RI-SMOM scheme. We will study the local current constructed from HISQ quarks since this is the simplest current and used in many analyses, such as the connected hadronic vacuum polarisation contribution to the anomalous magnetic moment of the muon [33].
In [11] the renormalisation constant for the HISQ local current was calculated using the form factor method discussed in Section II. Results are given for very coarse, coarse and fine lattices in Table IV of that reference. The calculation was done using valence s quarks and the form factor was determined for the local temporal vector current between two ss pseudoscalar mesons at rest 1 . From the discussion in Section II we expect such a determination of Z V to be exact so that Z loc V (F(0)) is equal to a perturbative series in α s that matches the lattice scheme to the MS scheme, up to discretisation effects. This was tested in [11] (Appendix B) by fitting the Z V results to this form, including the known O(α s ) coefficient in the perturbative series. A good fit was obtained that allowed values for Z loc V (F(0)) to be inferred on finer lattices. Here we will calculate Z loc V (SMOM) and compare it to Z loc V (F(0)). They should both contain the same perturbative series (since this is unique for a given operator) and differ only by discretisation effects.
To calculate Z loc V (SMOM) a little care is required in the construction of the SMOM vector vertex function with HISQ quarks. The operator / qq µ Λ µ V of Eq. (12) must be constructed to be a taste singlet. For a local (in spintaste notation, γ µ ⊗ γ µ ) current Λ µ V will have taste γ µ . GeV, plotted as a difference to the corresponding ZV at that lattice spacing obtained from the vector form factor at zero momentum-transfer. The fit shown (see text) accounts for discretisation errors and condensate contributions, but no condensate contributions are seen and they are strongly constrained to zero by the fit. The lower plot shows the same difference but for a Z loc V (SMOM) derived from results at µ = 2 GeV and 3 GeV in such a way as to reduce discretisation effects (see Eq. (30)). The fit here is to a simple (constant + a 4 ) form as described in the text.
This means that the / q in the vertex function must also have this taste. The correct construction is as Taking the spin-colour trace of this operator and dividing by 48q 2 then gives Z q /Z V . The wavefunction renormalisation is calculated in the same way as for the conserved current, Eq. (20). Results for Z loc V (SMOM) are given in Table III (column 3). This is after extrapolation to zero valence quark mass. Figure 5 shows that the impact of this is very small (we expect in this case that the mass dependence is purely a discretisation effect). Figure 6 (top plot) shows our results as a difference between Z loc V (SMOM) and Z loc V (F(0)). Z loc V (F(0)) values are from [11] and obtained on the same gluon field configurations that we use here. We plot the difference for the multiple µ values used for the Z loc V (SMOM) determination as a function of lattice spacing in Fig. 6. Results for a variety of µ values (given in column 2) on gluon field configurations at different lattice spacing values (denoted by the set number in column 1). Column 3 gives results using the RI-SMOM scheme and column 4 gives results using the standard RI -MOM scheme. Note that the RI -MOM results include the additional renormalisation factor of Eq. (25) ( Table II) that is needed to take the lattice current all the way to the MS scheme. Results are extrapolated to zero valence quark mass. Columns 5 and 6 give results for the modified (denoted by Rc) RI -MOM and RI-SMOMγ µ schemes in which a ratio to the value for the conserved current renormalisation in that scheme has been taken (Eq. (32)).  (30) are shown from very coarse to superfine lattice spacings noting that higher µ values are only accessible on finer lattices because of the constraint in Eq. (17). We can readily fit this difference of Z loc V values, ∆Z loc V , to a function constructed from possible discretisation effects. To keep the fit as general as possible we also allow for the existence of condensate terms to see to what extent they are constrained by the fit. We also allow for condensate terms multiplied by discretisation effects that would vanish in the continuum limit (and are therefore benign). We use cond,a 2 (aΛ/π) 2 ] .
All coefficients are given priors 0 ± 1. This fit has a χ 2 /dof value of 0.18 and finds no significant condensate contribution. The lowest order (1/µ 2 ) condensate term is constrained by the fit to have a very small coefficient compatible with zero: -0.020(44) (compare. Eq. (27)).
Thus we see that ∆Z loc V is compatible with being, as expected, purely a discretisation effect.
We have shown here that the Z V value obtained for the nonconserved local HISQ current using the RI-SMOM scheme is indeed exact i.e. it has no nonperturbative condensate contributions (visible at our high level of accuracy) that would survive the continuum limit as a source of systematic error. This can be traced to the fact that the condensate contributions present in the vector vertex function for the conserved vector current and in the inverse propagator must cancel because of the Ward-Takahashi identity. This identity also protects Z V from any effects arising from the gauge-fixing procedure.
This means that there is in fact no lower limit in principle to the µ value that can be used for the vector current renormalisation in the RI-SMOM scheme. In Fig. 6 (top plot) we include values corresponding to µ = 1 GeV. These show smaller discretisation effects than those for the higher µ values and so may be preferable on these grounds if only one µ value is used (which is all that is necessary in principle since no allowance needs to be made for condensate effects). The statistical errors possible with 20 configurations grow as µ is reduced. However, for µ = 1 GeV the uncertainties could still readily be reduced to the 0.1% level with higher statistics.
Smaller discretisation effects are possible by extrapolating in µ to µ = 0. A simple method that removes µ 2 a 2 terms in Z loc V (SMOM) combines results at two different µ values (for a given lattice spacing) to determine a new value This can always be done, given that Z loc V (SMOM) only depends on µ through discretisation effects. We use µ 1 = 3 GeV and µ 2 = 2 GeV and Eq. (30) returns a precise result because the statistical uncertainties are very small on these µ values. We show the results of taking a difference to Z loc V (F(0)) for this new Z V value in the lower plot of Fig. 6. The points clearly have smaller discretisation effects compared to the original Z V values that they were derived from. Given that the discretisation effects in Z loc V (F(0)) were relatively small [11] we interpret this as a reduction of discretisation effects in Z loc V (SMOM). We can fit the points in the lower plot of Fig. 6 to a very GeV, plotted as a difference to the corresponding ZV at that lattice spacing obtained from the vector form factor at zero momentum-transfer. The fit shown (see text) accounts for discretisation errors and condensate contributions, with condensate contributions being necessary to obtain a good fit.
simple curve : C + D(a × 1GeV) 4 and C is found to be 0.00008(66). The smaller discretisation effects seen using Eq. (30) may make this approach preferable to that of using Z loc V (SMOM) for a single µ value although it doubles the cost. Using three values of µ a higher-order scheme could obviously be devised to reduce discretisation effects further.

E. ZV for the local current in the RI -MOM scheme
We now turn to the determination of the renormalisation constant for the nonconserved local vector current using the RI -MOM scheme, Z loc V (MOM). Again, the vector vertex function must be a taste-singlet. The RI -MOM scheme uses a simple γ µ projector (Eq. (16)), which for the HISQ local vector current needs to have spin-taste γ µ ⊗ γ µ . Then we use to determine Z V /Z q along with Eq. (20) to determine Z q . Figure 7 shows the valence mass extrapolation for one set of raw results. Despite having a more significant mass extrapolation than for the RI-SMOM results (Fig. 5), this is still very mild. Table III gives our results in column 4, where we note that the values given for Z loc V (MOM) include the additional renormalisation factor shown in Eq. (25) and given in Table II. Figure 8 shows our results given, following the discussion in Section IV D, as a difference to the renormalisation constants obtained for the local current using the form factor method in [11]. This figure is very different from Fig. 6, with the results showing no sign of converging to zero in the continuum limit that would demonstrate agreement between the form factor and RI -MOM schemes for Z V . This shows the presence of condensate contributions in Z loc V (MOM) and to fit these results we need to include condensates that survive the continuum limit in the fit form.
For the difference of Z loc V values shown in Fig. 8 we use the same fit form as that used earlier for the RI-SMOM results in Eq. (29) (with the addition of an α 4 s to allow for uncertainty in the matching from MOM to MS as used in Eq. (26); this term has very little effect). This fit, with χ 2 /dof =0.14 is shown by the dashed lines in Fig. 8. It returns a coefficient for the leading-order condensate term of -0.209(63) which is consistent with the leading-order condensate term seen in the conserved current Z V calculated in the RI -MOM scheme (Eq. (27), with opposite sign because of our definition of ∆Z V here). Note the difference with the results in the RI-SMOM case.
The results of Fig. 8 show that the standard RI -MOM scheme cannot be used to determine an accurate result for Z V for nonconserved currents. If no attention is paid to the contamination of Z V by condensate contributions then O(1%) systematic errors will be made.
We can modify the RI -MOM scheme to address this issue, however. We know that the conserved current and the renormalised local current are the same operator in the continuum limit and so their vertex functions must contain the same nonperturbative contributions from the RI -MOM scheme in that limit. We can therefore calculate Z loc V (MOM) by taking a ratio of the vertex functions of the local and conserved currents. We call this scheme the RI -MOM Rc scheme. Specifically we calculate .
Taking the ratio also means that no additional renormalisation is needed in this case.
Our results from implementing this scheme are given in Table III (column 5). Figure 9 shows the results given . Z loc V (MOMRc) from the modified RI -MOM scheme of Eq. 32, plotted as a difference to the corresponding ZV at that lattice spacing obtained from the vector form factor at zero momentum-transfer. Results are shown for µ values from 2 GeV to 4 GeV. The fit shown (see text) accounts for discretisation errors and condensate contributions, but condensate contributions are strongly constrained to be zero.
once again as a difference to the renormalisation constant obtained for the local current in the form factor method. We now see that the difference of Z V values clearly approaches 0 in the continuum limit and there is no sign of condensate contamination in that limit. The results in the RI -MOM Rc scheme look very similar to those in the RI-SMOM scheme (see Fig. 6). We can fit the values for ∆Z loc V in Fig. 9 to the same form as that used for the RI-SMOM results (Eq. (29)). The fit gives χ 2 /dof = 0.32 and constrains the lowest-order condensate coefficient that would survive the continuum limit to -0.01 (5).
We conclude that the modified RI -MOM scheme of Eq. (32) does provide a method to determine an accurate renormalisation for the local vector current. The method does require calculations with the conserved current and so is more complicated than the RI-SMOM scheme.
F. ZV for the local current in the RI-SMOMγ µ scheme An alternative momentum-subtraction scheme is the RI-SMOM γµ scheme which uses the same vertex function (and wavefunction renormalisation) as the RI -MOM scheme but uses RI-SMOM kinematics (i.e. q = p 1 −p 2 = 0, p 2 1 = p 2 2 = q 2 = µ 2 ). To obtain an accurate result for Z V for the local current (as an example of a nonconserved current) we must modify the scheme as was done for the RI -MOM scheme in Eq. (32). The only difference is that we must also modify the tree-level vertex function factor for the conserved current from that of Eq. (24) to reflect the SMOM kinematics. Table III gives our results from this modified RI-SMOM γµ,Rc scheme in column 6. Figure 10 plots the difference of these Z V values with those from using the form factor method. We see that, as for the SMOM V (SMOMγ µ,Rc ) from the modified RI-SMOMγ µ scheme, plotted as a difference to the corresponding ZV at that lattice spacing obtained from the vector form factor at zero momentum-transfer. Results are shown for µ values from 2 GeV to 4 GeV. The fit shown (see text) accounts for discretisation errors and condensate contributions, but condensate contributions are strongly constrained to be zero. scheme in Fig. 6 and the modified RI -MOM scheme in Fig. 9, the values converge to zero as a → 0 as discretisation effects should. Discretisation effects are significantly larger here than in the previous schemes, however. We fit the results to the same functional form as used for the other schemes (i.e. Eq. (29)) and obtain a good fit (we double the prior width on (aµ) n terms to allow for the larger discretisation effects). χ 2 /dof = 0.56 and the lowest order condensate coefficient is constrained very tightly, as in the other exact cases, to -0.03 (5).
The same conclusions apply as for the RI -MOM scheme, i.e. that defining Z V from the ratio of vertex functions with the conserved current gives an exact result.

G. Renormalisation of the axial vector current
The renormalisation factors for axial vector currents can also be calculated using momentum-subtraction schemes. However, for actions with sufficient chiral symmetry the axial vector current renormalisation, Z A , can be related to the vector current renormalisation at zero quark mass. For example, for staggered quarks, Z S⊗T = Z S5⊗T 5 to all orders in perturbation theory [34]. Here S ⊗ T indicates the operator spin-taste and S5 = γ 5 S. This means that the local axial vector current and local vector current have the same renormalisation factor.
Having shown that the local vector current renormalisation factor can be calculated accurately and without contamination by condensate contributions in the RI-SMOM scheme in Section IV D, it therefore makes sense to use this value also for the local axial vector current. Indeed, doing a separate calculation of Z loc A risks introducing condensate contributions where none would be found using Z A = Z V . Figure 11 shows the difference  11. The difference between the vertex functions for the local vector and local axial vector currents as a function of µ in the RI -MOM and RI-SMOM schemes. The RI-SMOM scheme gives a much smaller nonperturbative contribution to ZA −ZV than the RI -MOM scheme, but neither scheme gives ZV = ZA which is true to all orders in perturbation theory. The values plotted result from an extrapolation to zero valence quark mass and are shown for the finest lattice we use (set 7 in Table I).
between the local vector and local axial vector vertex functions after extrapolation to zero quark mass, on the superfine lattices, set 7. Each point plotted is the difference of the local vector and local axial vector vertex We see that the difference in the RI-SMOM scheme is small but not zero. The results demonstrate approximately µ −6 behaviour expected on the basis of a chiral symmetry breaking condensate contribution [14]. Note that this contribution comes from Z q /Z A . For the RI -MOM scheme the difference is much larger then for RI-SMOM and has a smaller slope in this log-log plot. This reflects the known impact of chiral symmetry breaking nonperturbative artefacts in this scheme [14]. In both cases it would be preferable to use Z A = Z V , in the RI -MOM case using the modified RI -MOM Rc approach of Eq. (32).

V. INCLUDING QUENCHED QED EFFECTS
As lattice QCD calculations reach sub-percent precision it will become necessary to evaluate the electromagnetic corrections expected at this level. If QED effects are included in calculations involving nonconserved vector currents, such as the ongoing Fermilab/HPQCD/MILC calculations of the hadronic vacuum polarisation contribution to the anomalous magnetic moment of the muon [35], then consistency requires that QED effects are also included in the vector current renormalisation. Here we will study the impact of the valence quarks having electric charge on the renormalisation of the local vector current using the RI-SMOM scheme (for earlier results using different methods see [36,37] FIG. 12. Results for the renormalisation factor, ZV for the QCD+QED conserved current for the HISQ action, calculated using the RI-SMOM scheme. Results are given for coarse, fine and superfine gluon field configurations for quark electric charge, Q = 2e/3 and a variety of momenta with magnitude aµ in lattice units. We include 'quenched QED' in our lattice calculations by multiplying our QCD gauge fields by a U(1) gauge field representing the photon. The photon field, A µ (k), is randomly generated in momentum space from a Gaussian distribution with variance 1/k 2 to yield the correct O(a 2 )-improved Feynman gauge propagator on the lattice (the definition ofk is given in Eq. (21)). A µ (k) is then converted to Landau gauge and transformed to position space. To make sure of the correct gauge covariance in position space it is important to remember that the position of the gauge fields is at the centre of the links, and not the sites [38]. The A µ field in position space is then used as the phase to construct a U(1) field [39] in the form exp(ieQA µ ) where Q is the charge of the quark that will interact with the field, in units of the charge on the proton, e. We use the QED L formulation of compact QED [40], in which all zero modes are set to zero, A µ (k 0 , k = 0) = 0 with A µ in Landau gauge (for a review of approaches to handling zero modes in QED on the lattice see [41]). We multiply the gluon field for each link of the lattice by the appropriate U(1) field before applying the HISQ smearing. The valence quarks can then interact with the photon via the standard HISQ action. Note that the sea quarks remain electrically neutral, so this is not a fully realistic scenario. Nevertheless it allows us to evaluate the most important QED effects.
We have tested that the U(1) configurations we generate correctly reproduce the O(α QED ) perturbation theory prediction for the average plaquette [42], independent of gauge choice. Our results for the average value of the U(1) link field also agree with the O(α QED ) expectations: Landau gauge : 1 − α QED Q 2 0.0581 (33) Feynman gauge : 1 − α QED Q 2 0.0775 Note that the Landau gauge O(α QED Q 2 ) coefficient is 1/C F = 3/4 that of the corresponding QCD result for the a 2 -improved gluon action [43] since the gluon propagator then has the same form as that of the photon here. The Feynman gauge coefficient is then 4/3 of the Landau gauge coefficient.
Although we have tested calculations as a function of quark charge, Q, the results we will show here are all for Q = 2/3. The results are not extrapolated to zero valence quark mass and are instead just the values at the sea light quark mass on each ensemble. The valence mass dependence of the results is observed to be negligibly small, as was the case in pure QCD.
An important test of the interaction between the quarks and the QCD+QED gauge fields is that Z V = 1 for the QCD+QED conserved current in the RI-SMOM scheme, as expected from a trivial extension of the Ward-Takahasi Identity to this case. This is demonstrated in Fig. 12.
Our analysis for the renormalisation of the local vector current in the RI-SMOM scheme will study the ratio of the Z loc V calculated with and without the inclusion of electromagnetic effects. We proceed exactly as for the pure QCD case discussed in Section IV D. The strong correlations between the QCD and QCD+QED calculations allow very precise determination of this ratio (a typical correlation being ∼ 0.99). We will denote a quantity X calculated in pure QCD as X[QCD] while the same quantity calculated with the inclusion of QED effects will be denoted X[QED + QCD]. We will also employ the notation Because QED is a long-range interaction it is important to test finite-volume effects, although we do not expect them to be large here since we studying renormalisation of electrically neutral currents. The finite-volume effects in the self-energy function of fermions has been studied in [44] with the result that for off-shell quarks the finite-volume effects start at order 1/L 2 s where L s is lattice spatial extent. This implies that even the finitevolume effects for quantities such as Z q should be small. Figure 13 confirms both of these expectations with results on the three lattice sets with the same parameters but different volumes (sets 3, 4 and 5, ranging in spatial extent from 2.9 fm to 4.9 fm). Negligible effects are seen here and we therefore ignore finite-volume issues in the following analysis.
Our results for the effect of quenched QED on Z V for the local HISQ current in the RI-SMOM scheme are given for µ values from 2 GeV to 4 GeV and at 3 values of the lattice spacing in Table IV. The results are plotted in Fig. 14.
Given our results for the pure QCD case in Section IV D we expect the results for Z V for QCD+QED to be similarly well-behaved. We therefore perform a fit to the ratio of Z V for QCD+QED to that for pure QCD results that allows for both discretisation effects along with a perturbative expansion for the ratio of renormalisation constants. The leading QCD effects will cancel between the numerator and denominator of the ratio and so the leading term in this expansion will be O(α QED ).  We can even fix the coefficient of the leading-order term based on the QCD perturbation theory for the pure QCD case. The O(α s ) coefficient for Z loc V for pure QCD is -0.1164(3) [11]. We therefore expect that the coefficient of α QED Q 2 in the QED case is −0.1164 × 3/4 = −0.0873. Z loc V of 0.9997, very close to 1. There will be in principle α s α QED corrections to this which are likely to have an even smaller impact.
We therefore take a fit form for the ratio of Z V values given in Table IV  We use i = 0, 1, 2, 3 and j = 1, 2, 3 fixing c 0 to the value given above. Note that α QED does not run in this expression because we are using quenched QED. α s in Eq. (34) is taken as α MS (1/a). This fit returns a χ 2 /dof value of 0.25. The fit is plotted with the results in Fig. 14. Figure 14 shows that the results for Z V behave as expected. The impact of quenched QED on the value of Z loc V is tiny and indeed negligible if we imagine working to an accuracy of 0.1%. Note that this follows directly from the analysis above in which we derive the O(α QED ) coefficient for the QCD+QED case from the pure QCD case. Because the HISQ action is so highly improved Z loc V is very close to 1 in the pure QCD case. It then has to be true that the difference from 1 in Z V induced by QED will be over 100 times smaller than that induced by QCD. For the HISQ action this means that the impact of QED in Z loc V is of order 0.03%. This should be contrasted to the case from the domain-wall action where the Z V value in pure QCD is 0.7 and so the impact of quenched QED is to change Z V by approximately 0.3/100 for Q = 2e/3, in this case 0.3% (see Table 6 of [36]); this is not negligible.
The effect of having electrically charged sea quarks would appear in Z V at O(α 2 s α QED ) i.e. two orders in α s below the leading term; the leading term comes from a photon exchange across a quark bubble created from a gluon. This is unlikely to change the picture significantly. The effect of QED on Z V is of course not a physical result and it needs to be combined with hadronic matrix elements for the vector current to understand the physical effect of QED. For this we simply take the values for Z V at a fixed µ value for the ensembles for which we have matrix element results, multiply them and extrapolate to the continuum limit. Different quark formalisms should agree on the physical effect (on an uncharged sea). We will give an analysis of the impact of quenched QED on vector current matrix elements calculated with the HISQ action elsewhere.

VI. CONCLUSIONS
We have shown by explicit calculation how the vector Ward-Takahashi identity works for the HISQ action in lattice QCD. Renormalisation methods that make use of this identity will give a renormalisation constant of 1 for the conserved current as would be obtained in continuum QCD. The RI-SMOM momentum-subtraction scheme is such a scheme but the RI -MOM scheme is not and this has implications for the accuracy achievable for Z V for nonconserved currents within each scheme. Our calculations have used the HISQ action but our conclusions are not specific to this action.
The RI-SMOM scheme provides precise values for Z V for nonconserved currents (using momentum-sources) that are completely nonperturbative. Our results show that the Z V values are 'exact' in being free of condensate contamination. This means that we can simply determine Z V at a given momentum scale µ on a given gluonfield ensemble, multiply our vector current hadronic matrix element by it and then extrapolate results for the renormalised matrix element to the continuum limit. Because there is no condensate contamination there is no lower limit to the µ value that can be used. Statistical errors grow as µ is reduced but discretisation effects become smaller. In Section IV D we demonstrated a simple method to reduce discretisation effects, if they are an issue, by combining results from two different µ values.
The RI -MOM scheme can also provide precise values for Z V for nonconserved currents, but is not completely nonperturbative. A more critical problem with this scheme is that the Z V values for both conserved and nonconserved currents have condensate contributions that begin at 1/µ 2 . This means that the Z V values cannot be used to obtain accurate renormalised vector current matrix elements in the continuum limit without an analysis of these condensate contributions. This requires numbers for Z V at multiple µ values and a fit that includes condensate terms. If this analysis is not done, the results obtained in the continuum limit will be incorrect at the 1% level.
An alternative to the standard RI -MOM scheme that avoids this problem is to determine Z V from a ratio of vector vertex functions for the conserved and nonconserved currents. We call this scheme RI -MOM Rc . A similarly modified RI-SMOM γµ scheme can also be used to obtain an exact Z V . These schemes are discussed in Sections IV E and IV F. It is straightforward to include quenched QED effects in the determination of the vector current renormalisation factor in a fully nonperturbative way using the RI-SMOM scheme and to obtain a full understanding of the results (including consistency with perturbation theory). We see only very small (below 0.1%) effects for the local HISQ vector current reflecting the fact that the renormalisation factors in the pure QCD case are already very close to 1. We will include the QCD+QED Z V values in a future QCD+QED determination of hadronic vector current matrix elements.  B1)). This is the mean value of the gluon field Uµ in Landau gauge. Column 3 gives the results for the ZV values determined from the form factor using the matrix element of the temporal 1link current between two pions at rest a . The asterisk next to set 6 is to denote that the results given here are actually for another fine ensemble with am l = 0.0074, ams = 0.037 and amc = 0.44.
Koponen and A.C. Zimermmane-Santos for providing the Z V (F(0)) values. in conjunction with higher-order difference operators for ∆ µ,± but we do not do that here.

Appendix B: Renormalisation of the 1-link vector current
Quark-line disconnected contributions for vector current-current correlators require the use of a tastesinglet vector current for staggered quarks. This has the same taste as the conserved current but it is often more convenient to use a simpler current than the conserved one. Here we discuss the renormalisation of the nonconserved 1-link point-split vector current using momentumsubtraction schemes. The qualitative picture is the same as that for the local current and so we simply include RI-SMOM results in this Appendix for completeness. They are relevant to our ongoing calculations of, for example, the quark-line disconnected pieces of the hadronic vacuum polarisation contribution to the anomalous magnetic moment of the muon. We consider the 1-link point-split vector current with spin-taste (γ µ ⊗ I). The operators that we use include gluon fields between the point-split quark fields to maintain gauge invariance. We take these gluon fields to be 'thin links' i.e. no smearing is applied to them. We considered the two simplest constructions of this current. One, which we denote the forward 1-link operator, is the conserved current with all 3-link terms removed: The other 1-link operator we consider is the symmetric operator j 1link-symm µ ≡ 1 2u 0 ψ(x)(γ µ ⊗ 1)U µ (x)ψ(x +μ) (B2) The two definitions coincide with the MOM kinematics. In the SMOM case, while the values produced from the two different definitions are not identical they agree within our statistical uncertainties. In what follows we then only present results for the forward 1-link current. Note that in the definitions of the 1-link current above we include a factor 1/u 0 . u 0 is a 'tadpole-improvement' factor [45] which can be used, as here, to reduce the mismatch between lattice currents containing gluon fields and their continuum counterparts. u 0 works by cancelling universal effects from tadpole diagrams that arise from the construction of the lattice gluon field. u 0 can in principle be any suitable ensemble average of a function of the gluon field that achieves this. Here we use the mean value of the gluon field U µ in Landau gauge as the most appropriate form of u 0 in this case. The values for u 0 depend on the ensemble and are listed in Table V 2 . 2 Z V for the tadpole-improved current is u 0 times Z V for the cur-We proceed for the 1-link case in the same way as for the cases discussed in the main body of the paper. The wavefunction renormalisation is exactly the same as before. We calculate the vertex function for the 1-link current using an appropriate projector. For the RI-SMOM case we use aq ν (aq) 2 Tr (γ ν ⊗ I)Λ µ V .
(B3) In determining Z 1link V an additional technical detail for point-split operators is that we must divide the vertex function in the full theory by the result of the treelevel (noninteracting) case. This was discussed previously for the conserved current in the RI -MOM case in Section IV C (and denoted V γ⊗I in Eq. (23) and above). The tree-level result for the forward 1-link current for the RI-SMOM kinematics is: γ⊗I (SMOM) = 1 2 µ (e iap2,µ(S−T )µ + e −iap1,µ(S−T )µ ) (B4) A further technical detail arises when using twisted boundary conditions to insert momentum with pointsplit operators in the vertex functions that we calculate. The propagator with twisted momentum can be written in terms of the untwisted one as S(x, p) = e −iθx S(x, p + θ). (B5) We want the vertex function for a point-split operator to take the following form (using a 1-link operator Γ as an example, but dropping the gluon fields for clarity): x γ 5 e i(p1+θ1)x S † (x, p 1 + θ 1 ) (B6) × γ 5 Γ µ e −i(p2+θ2)x S(x +μ, p 2 + θ 2 ) = x γ 5 e ip1xS † (x, p 1 )γ 5 Γ µ e −ip2x e iaθ2,µS (x +μ, p 2 ).
The factor e iaθ2,µ has to be inserted by hand.
Our results for Z 1link V (SMOM) using the RI-SMOM scheme are given in Table VI for a variety of µ values for three values of the lattice spacing. We expect the Z V values obtained with RI-SMOM to be well-behaved and free of condensate contributions because of the protection of the Ward-Takahashi identity, as for the local current discussed in Section IV D. We can test this, as was done for the local case, by taking a difference of the Z V values with those obtained from the form factor method.
The results for Z V from the form factor method are given in Table V for a variety of µ values and on ensembles with a range of lattice spacing values. The results for the difference of Z V values between the form factor and rent with no tadpole-improvement.
RI-SMOM methods is plotted in Fig. 15. We show the results of a simple fit to a sum of possible discretisation effects: Here α s is in the MS scheme at scale 1/a. We have to include (aΛ) n terms as well as (aµ) n terms here because of the relatively large discretisation effects in the Z V values obtained from the form factor method. The priors on the coefficients of the fit are taken as: 0 ± 3. The fit gives a χ 2 /dof = 0.9. This confirms that, again in this case, the RI-SMOM method gives a well-behaved result for Z V .
The conclusion from this is that the renormalisation factors for the 1-link current obtained in the RI-SMOM scheme on the lattice can be used straightforwardly, and in a fully nonperturbative way, to renormalise matrix elements of the 1-link current obtained in a lattice calculation. This means that values can be taken, for example from Table VI, for a fixed µ value on each ensemble. The µ chosen can take any value and the only limitation on taking it to have a small value (for minimal discretisation effects) is that of the statistical errors that grow as µ is reduced.