Strange and charm HVP contributions to the muon (g − 2) including QED corrections with twisted-mass fermions

We present a lattice calculation of the Hadronic Vacuum Polarization (HVP) contribution of the strange and charm quarks to the anomalous magnetic moment of the muon including leading-order electromagnetic corrections. We employ the gauge configurations generated by the European Twisted Mass Collaboration (ETMC) with Nf = 2 + 1 + 1 dynamical quarks at three values of the lattice spacing (a ≃ 0.062, 0.082, 0.089 fm) with pion masses in the range Mπ ≃ 210-450 MeV. The strange and charm quark masses are tuned at their physical values. Neglecting disconnected diagrams and after the extrapolations to the physical pion mass and to the continuum limit we obtain: aμs(αem2) = (53.1 ± 2.5) · 10− 10, aμs(αem3) = (−0.018 ± 0.011) · 10− 10 and aμc(αem2) = (14.75 ± 0.56) · 10− 10, aμc(αem3) = (−0.030 ± 0.013) · 10− 10 for the strange and charm contributions, respectively.


Introduction
The anomalous magnetic moment of the muon a µ ≡ (g − 2)/2 is known experimentally with an accuracy of the order of 0.54 ppm [1], while the current precision of the Standard Model (SM) prediction is at the level of 0.4 ppm [2]. The tension of the experimental value with the SM prediction, a exp µ − a SM µ = (28.8 ± 8.0) · 10 −10 [2], corresponds to 3.5 standard deviations and might be an exciting indication of new physics.
The forthcoming g − 2 experiments at Fermilab (E989) [3] and J-PARC (E34) [4] aim at reducing the experimental uncertainty by a factor of four, down to 0.14 ppm. Such a precision makes the comparison of the experimental value of a µ with theoretical predictions one of the most important tests of the Standard Model in the quest for new physics effects.
It is clear that the experimental precision must be matched by a comparable theoretical accuracy. With a reduced experimental error, the uncertainty of the hadronic corrections will soon become the main limitation of this test of the SM. For this reason an intense research program is under way to improve the evaluation of the leading order hadronic contribution to a µ due to the Hadronic Vacuum Polarization (HVP) correction to the oneloop diagram, a had µ (α 2 em ), as well as to the next-to-leading-order hadronic ones. The latter include the O(α 3 em ) contribution of diagrams containing HVP insertions and the leading hadronic light-by-light (LBL) term [5].

JHEP10(2017)157
The theoretical predictions for the hadronic contributions have been traditionally obtained from experimental data using dispersion relations for relating the HVP function to the experimental cross section data for e + e − annihilation into hadrons [6,7]. An alternative approach was proposed in refs. [8][9][10], namely to compute a had µ (α 2 em ) in Euclidean lattice QCD from the correlation function of two electromagnetic currents. In this respect an impressive progress in the lattice determinations of a had µ (α 2 em ) has been achieved in the last few years [11][12][13][14][15][16][17][18][19][20][21] and very interesting attempts to compute also the LBL contribution are under way both on the lattice [22,23] and via dispersion approaches and Chiral Perturbation Theory (ChPT) [24][25][26].
With the increasing precision of the lattice calculations, it becomes necessary to include electromagnetic (e.m.) and strong isospin breaking (IB) corrections (which contribute at order O(α 3 em ) and O(α 2 em (m d − m u )), respectively) to the HVP. In this paper we present the results of a lattice calculation of the leading radiative e.m. corrections to the HVP contribution due to strange and charm quark intermediate states, obtained using the expansion method of refs. [27,28]. Given the large statistical fluctuations, we are not in the position of giving results for the e.m. and IB corrections to the HVP contribution from the light up and down quarks, although we will give some details of our computation. For the same reason we do not have yet results for the disconnected contributions.
The main results of the present study for a had Our findings demonstrate that the expansion method of refs. [27,28], which has been already applied successfully to the calculation of e.m. and strong IB corrections to meson masses [28,29] and leptonic decays of pions and kaons [30,31], works as well also in the case of the HVP contribution to a µ . This is reassuring about the feasibility of the determination of the leading e.m. and strong IB corrections to the HVP contribution from the up and down quarks, which is expected to be not negligible [5] and will be addressed in a separate work. For a recent calculation of these corrections, though at a large pion mass and at a fixed lattice spacing, see ref. [32], where, as expected, the strong IB effect is found to be at the percent level. In the strange and charm sectors the e.m. corrections (1.3)-(1.4) are found to be negligible with respect to present uncertainties. The paper is organized as follows. In section 2 we introduce the basic quantities and notation. In section 3 we describe the lattice calculation and give the simulation details. In section 4 we present the calculation of the strange and charm contributions to the HPV at order O(α 2 em ) and in section 5 the corresponding e.m. corrections at order O(α 3 em ), which represent the original part of this work. Finally, section 6 contains our conclusions and outlooks for future developments. can be related to the Euclidean space-time HVP function Π(Q 2 ) by [8][9][10] where Q is the Euclidean four-momentum and the kinematical kernel f (Q 2 ) is given by with m µ being the muon mass and ω ≡ Q/m µ . The HVP form factor Π(Q 2 ) is defined through the HVP tensor as where . . . means the average of the T -product of the two electromagnetic (e.m.) currents over gluon and fermion fields and with q f being the electric charge of the quark with flavor f in units of e. In eq. (2.1) the subtracted HVP function Π R (Q 2 ) ≡ Π(Q 2 ) − Π(0) appears. This is due to the fact that the photon wave function renormalization constant absorbs the value of the photon self energy at Q 2 = 0 in order to guarantee that the e.m. coupling α em is the experimental one in the limit Q 2 → 0.
The HVP function Π R (Q 2 ) can be determined from the vector current-current Euclidean correlator V (t) defined as Taking into account that V (−t) = V (t) and choosing Q along the time direction only, one has [33] Consequently the HVP contribution a had µ can be written as wheref (t) is given by [33]

JHEP10(2017)157
In what follows we will limit ourselves to the connected contributions to a had µ . In this case each quark flavor f contributes separately, i.e.
For sake of simplicity we drop the label f and the suffix (conn), but it is understood that hereafter we refer to the connected part of a had µ for a generic quark flavor f .
The vector correlator V (t) can be calculated on a lattice with volume L 3 and temporal extension T at discretized values of t ≡ t/a from −T /2 to T /2 with T = T /a. From now on all the "overlined" quantities are in lattice units. A natural procedure is to split eq. (2.7) into two contributions, a had µ (<) and a had µ (>), corresponding to 0 ≤ t ≤ T data and t > T data , respectively. In the first contribution a had µ (<) the vector correlator is directly given by the lattice data, while for the second contribution a had µ (>) an analytic representation is required (see refs. [17,18,20,21]). If T data is large enough that the ground-state contribution is dominant for t > T data , one can write where is the (squared) matrix element of the vector current operator (for the given quark flavor f ) between the vector ground-state and the vacuum. Note that T data cannot be taken equal to T /2, because on the lattice the vector correlator possesses backward signals. In order to avoid them one has to choose an upper limit T data sufficiently smaller than T /2. An important consistency check is that the sum of the two terms in the r.h.s. of eq. (3.1) should be almost independent of the specific choice of the value of T data , as it will be shown later in section 4.

Simulation details
The gauge ensembles used in this work are the same adopted in ref. [34] to determine the up, down, strange and charm quark masses. We employed the Iwasaki action [35] for gluons and the Wilson Twisted Mass Action [36][37][38] for sea quarks. In order to avoid the mixing of strange and charm quarks in the valence sector we adopted a non-unitary set  [34]). The values of the strange and charm quark bare masses aµ s and aµ c , given for each gauge ensemble, correspond to the physical strange and charm quark masses determined in ref. [34]. The central values and errors of the pion, kaon and D-meson masses are evaluated using the bootstrap events of the eight branches of the analysis of ref. [34].
up [39] in which the valence strange and charm quarks are regularized as Osterwalder-Seiler fermions [40], while the valence up and down quarks have the same action of the sea. Working at maximal twist such a setup guarantees an automatic O(a)-improvement [38,39]. We considered three values of the inverse bare lattice coupling β and different lattice volumes, as shown in table 1, where the number of configurations analyzed (N cf g ) corresponds to a separation of 20 trajectories. At each lattice spacing, different values of the light sea quark masses have been considered. The light valence and sea quark masses are always taken to be degenerate. The bare masses of both the strange (aµ s ) and the charm (aµ c ) valence quarks are obtained, at each β, using the physical strange and charm masses and the mass renormalization constant (RC) determined in ref. [34]. The values of the lattice spacing are: a = 0.0885(36), 0.0815(30), 0.0619(18) fm at β = 1.90, 1.95 and 2.10, respectively.
In this work we made use of the bootstrap samplings elaborated for the input parameters of the quark mass analysis of ref. [34]. There, eight branches of the analysis were adopted differing in: • the continuum extrapolation adopting for the scale parameter either the Sommer parameter r 0 or the mass of a fictitious PS meson made up of strange(charm)-like quarks; • the chiral extrapolation performed with fitting functions chosen to be either a polynomial expansion or a Chiral Perturbation Theory (ChPT) Ansatz in the light-quark mass; • the choice between the methods M1 and M2, which differ by O(a 2 ) effects, used to determine in the RI'-MOM scheme the mass RC Z m = 1/Z P .
The kernel function f (t), appearing in eqs. (3.2)-(3.3), is explicitly given by and can be easily calculated numerically at any value of t for given values of the muon mass in lattice units, m µ ≡ am µ . This is shown in figure 1 in the case of the muon at the three values of the lattice spacing of the ETMC ensembles of table 1 (left panel) and for various values of the lepton mass (right panel) ranging from the µ to the τ mass. The kernel function f (t) is proportional to t 4 at small values of t, diverges as t 2 at large values of the time distance and has some sensitivity to the value of the lattice spacing. Instead it changes significantly with the mass of the lepton enhancing the role of the large-time behaviour of the vector correlator in the case of light leptons.

Local versus conserved vector currents on the lattice
The vector correlator V (t) can be calculated using either the lattice conserved vector current J C µ (x) or the local one J µ (x). The latter needs to be renormalized and in our twistedmass setup the local vector current for each quark flavor f is given by where, being at maximal twist, the renormalization is multiplicative through the renormalization constant Z V .

JHEP10(2017)157
The variation of the lattice action with respect to a vector rotation α V (x) of the quark fields, i.e. ψ(x) → e iq f α V (x) ψ(x) and ψ(x) → ψ(x) e −iq f α V (x) (for any quark flavor f ), provides the relevant Ward-Takahashi identity for the conserved current J C µ expressed in terms of the backward lattice derivative. In our twisted-mass setup one has According to the vector Ward-Takahashi identity the polarization tensor J C µ (x)J C ν (y) is not transverse because of the contact term arising from the vector rotation of the conserved current J C ν (y), which generates the backward lattice derivative of the tadpole operator and is power divergent as 1/a 3 . Thus, in the case of two conserved currents the transverse HVP tensor is defined as where the tadpole operator is explicitly given by On the contrary, in the case of one conserved and one local currents there is no contact term because the vector rotation of the local current (3.5) is zero. One gets which is transverse only with respect to the µ index (i.e., ∂ b µ Π CL µν (x, y) = 0, where ∂ b µ is the backward lattice derivative).
In the case of two local currents the polarization tensor J µ (x)J ν (y) is not transverse. The mixing pattern of the product of two local currents with all possible operators with equal and lower dimensions has been investigated for the twisted-mass setup in ref. [41]. The outcome is that at maximal twist one has where Π µν (Q) is the transverse polarization tensor, S 5 and S 6 are the vacuum expectation values of the dimension-5 and -6 terms of the Symanzik expansion of the twisted-mass action, m is the (twisted) quark mass and the quantities Z 1 , Z m , Z L and Z T are mixing coefficients. In the r.h.s. of eq. (3.10) the second and third terms do not depend on Q, while the fourth and fifth terms are Q-dependent. The former ones can be eliminated by considering the subtracted form where we have considered that Π µν (0) = 0 in the infinite volume limit [33]. Choosing Q along the time direction only with µ = ν = i = 1, 2, 3 one has Using eqs. (2.3) and (2.5) and taking into account that In the renormalized HVP function [Π(Q 2 ) − Π(0)] the term (Z L + Z T ) cancels out, so that eq. (2.6) is recovered and the O(a)-improvement of the renormalized HVP function is guaranteed.
In this work the local version of the vector current is adopted (see later eq. (3.23) in section 3.4).

Perturbative QCD (pQCD) and the behavior of V (t) at small t
The HVP function Π R (Q 2 ) obeys the (once subtracted) dispersion relation where R had (s) is related to the (one photon) e + e − annihilation cross section into hadrons, σ had (s), by with s being the center-of-mass energy and s thr = 4M 2 π . The pQCD prediction for R had (s) is known up to three loops including mass corrections [42]. Here we limit ourselves to the lowest order prediction, which, for each quark flavor f , reads as where m is the on-shell quark mass. Inserting eq. (3.16) into eq. (3.14) one obtains which exhibits a logarithmic divergence.

JHEP10(2017)157
In the continuum the vector correlator (2.5) can be obtained simply by taking the Fourier transform of the spatial components of the HVP tensor (2.3). Choosing Q along the time direction only, one gets Using eq. (3.14) one has Consequently, using the pQCD result (3.16) for R had (s) (including the quark mass threshold s thr = 4m 2 ) the pQCD prediction for V (t) is given by and therefore at small values of t the vector correlator V pQCD (t) is dominated by a mass-independent term, namely which represents also the vector correlator V pQCD (t) in the massless limit.
In figure 2 we compare the pQCD predictions (3.21) and (3.22) with the vector correlator V (t) obtained using the ETMC ensembles A30.32, B25.32 and D20.48, which share an approximate common value of the light-quark mass m 12 MeV and differ only in the values of the lattice spacing. It can be clearly seen that at small values of t the lattice data match nicely the (lowest order) pQCD prediction. The inclusion of the radiative corrections from ref. [42] leads to an effect of the order of ≈ 10%, which does not modify the quality of the agreement shown in figure 2.
A closer look to figure 2 shows that the matching with pQCD is present up to time distances of ≈ 1 fm (the agreement can be extended in the case of the strange vector correlator by including the corrections due to the strange quark mass), which corresponds to 1/Λ QCD with Λ QCD ≈ 300 MeV, i.e., the agreement is observed down to energy scales of the order of Λ QCD . One would expect that pQCD works at t 1/Λ QCD or Q Λ QCD . The fact that instead the matching appears to work at larger time distances is a nice manifestation of the quark-hadron dualityà la SVZ, which states that the sum of the contributions of the excited states is dual to the pQCD behaviour. Finally, it is interesting to estimate the contribution to a had µ coming from values of Q 2 larger than Q 2 max 1/a 2 , which for our lattice setup is always larger than 4 GeV 2 . Using the pQCD prediction (3.17) for the large Q 2 -behavior of Π R (Q 2 ), one gets: a had µ (Q 2 > 4 GeV 2 ) 1.3, 0.11, 0.06 (in units of 10 −10 ) in the case of the light, strange and charm contributions, respectively. The above findings represent only a small fraction of the uncertainties of the present lattice estimates of the three contributions to a had µ (see refs. [13][14][15]19]).
Alternatively we can check the change induced in the kernel function f (t) by cutting the upper integration limit in eq. (3.4) to ω max = Q max /m µ 1/(am µ ). Since in our lattice setup ω max 20, the kernel function f (t) changes at most by one part over 10 6 at small t in the case of the muon.

Ground-state identification
Our numerical simulations of the vector correlator V (t) have been carried out in the context of a more general project aiming at the determination of the e.m. and strong IB corrections to pseudoscalar meson masses and leptonic decay constants [43]. In this project the bilinear operators were constructed adopting opposite values of the Wilson r-parameter. Thus, instead of eq. (3.5) the evaluation of the vector correlator has been carried out using the following local current:   [34]), obtained using the local currents (3.5) (same r-parameters) and (3.23) (opposite r-parameters) for the ETMC gauge ensembles A30.32, B25.32 and D20.48, which share an approximate common value of the light-quark mass and differ by the values of the lattice spacing. The empty markers represent the value extrapolated in the continuum limit assuming a linear behavior in a 2 . figure 3, where the contribution of the strange quark to a had µ , evaluated using either the current (3.23) or the connected insertion of eq. (3.5), is shown as a function of a 2 for the three ensembles A30.32, B25.32 and D20.48, which share an approximate common value of the light-quark mass. It can be seen that the same continuum limit is reached using either currents, confirming that the difference is due to discretization effects of order O(a 2 ). Moreover, the absence of disconnected insertions in the current (3.23) implies that the "purely connected" vector correlator based on the current (3.5) is a well defined quantity and admits the hadron decomposition necessary for having the representation (3.3) (see also refs. [17,18] and therein quoted).
The statistical accuracy of the meson correlators is based on the use of the so-called "one-end" stochastic method [44], which includes spatial stochastic sources at a single time slice chosen randomly. Four stochastic sources (diagonal in the spin variable and dense in the color one) were adopted at first per each gauge configuration.
In the case of the light-quark contribution the signal-to-noise ratio does not allow to determine the ground-state mass M V and the corresponding matrix element Z V from the behavior of the vector correlator at large time distances. This is at variance with the case of the strange and charm contributions, as it is illustrated in figure 4, where it is also shown that discretization effects are sub-leading.     Thus, the identification of the ground-state is presently possible only in the case ofss andcc vector mesons. To improve the statistics we took a significative advantage by using the DD − αAMG solver [45], which has allowed us to increase by a factor of 5 the number of stochastic sources in the case of the strange quark. In this way we find that the quality of the plateaux, shown in figure 5, is acceptable in the strange sector and nice in the charm one. In the case of the light-quark contribution an increase of the statistics by a factor ≈ 20 is expected to be needed.
For each gauge ensemble the masses M V and the matrix elements Z V are extracted from a single exponential fit (including the proper backward signal) in the range t min ≤ t ≤ t max . The values chosen for t min and t max are collected in table 2. The results for the strange contribution to a had µ (<), a had µ (>) and their sum a had µ obtained adopting four choices of T data , namely: T data = (t min + 2), (t min + t max )/2, (t max − 2) and (T /2 − 4), are collected in table 3 for illustrative purposes in the case of few ETMC gauge ensembles.
The separation between a had µ (<) and a had µ (>) depends on the specific value of T data , as it should be, but their sum a had µ is almost independent of the choice of the value of T data in the range between t min and t max . This is also reassuring of the fact that the value of a had µ is not contaminated significantly by the presence of backward signals in the correlator V (t).
In the case of the charm contribution the value of a had µ (>) is always several orders of magnitude smaller than a had µ (<) and the latter turns out to be the same for all the four choices of T data .
Note that for T data = T /2 − 4 the contribution a had µ (>), which depends on the analytic representation In what follows all the four choices of T data will be employed in the various branches of our bootstrap analysis. The corresponding systematics is largely sub-dominant with respect to the other sources of uncertainties and it will not be given separately in the error budget.
The results obtained for the strange and charm contributions to a had µ are shown by the empty markers in figure 6. We observe a mild dependence on the light-quark mass, being driven only by sea quarks, and also small residual FSEs visible only in the case of the a had µ (<) 38.03 (28) a had µ (<) 40.83 (14) 43.18 (17) Table 3. Results for the strange contribution to a had µ (<), a had µ (>) and their sum a had µ , in units of 10 −10 , obtained assuming T data = (t min + 2), (t min + t max )/2, (t max − 2) and (T /2 − 4) for the ETMC gauge ensembles A40.24, A30.32, B25.32 and D15.48. Errors are statistical only. strange contribution. The errors of the data turn out to be dominated by the uncertainties of the scale setting, which are similar for all the gauge ensembles used in this work.
In ref. [13] a modification of the calculated a had µ at pion masses above the physical point has been proposed in order to weaken the pion mass dependence of the resulting a had µ for improving the reliability of the chiral extrapolation. Though the procedure of ref. [13] has been conceived mainly for the light contribution to a had µ , we have explored its usefulness also in the case of the strange and charm contributions. The proposal consists in multiplying the Euclidean 4-momentum transfer Q 2 by a factor equal to (M V /M phys V ) 2 in order to modify the Q 2 -dependence of the HVP function Π R (Q 2 ) without modifying its value at the physical point. One obtains the same effect in our master formulae by redefining the lepton mass as  [2] have been adopted for the physicalss andcc vector meson masses, respectively.
The expected advantage of the use of the effective lepton mass (4.1) comes from the fact that the kernel function, and therefore a had µ , depends only on the lepton mass in lattice units (see eq. (3.4)). Thanks to eq. (4.1), which will be referred to as the Effective Lepton Mass (ELM) procedure, the knowledge of the value of the lattice spacing is not required and therefore the resulting a had µ is not affected by the uncertainties of the scale setting. The drawback of the ELM procedure is instead represented by its potential sensitivity to the statistical fluctuations of the vector meson mass extracted from the lattice data.
The results obtained adopting the ELM procedure (4.1) in the case of the strange and charm contributions to a had µ are shown by the filled markers in figure 6, where the physical values for thess andcc vector masses have been taken from PDG [2] (namely, M (phys) V = 1.0195 and 3.0969 GeV, respectively). 1 It can be seen that the ELM procedure reduces remarkably the overall uncertainty of the data. Moreover, it further weakens the pion mass dependence (in any case driven only by the sea quarks) and modifies the discretization effects, leading to a better scaling behavior of the data in the case of the charm contribution. Since the pion mass dependence is in any case quite mild, the ELM procedure can be viewed as an alternative way to perform the continuum extrapolation and to avoid the scale setting uncertainties.
Using the data obtained either with or without the ELM procedure we have performed a combined fit for the extrapolation to the physical pion mass, the continuum and infinite 1 We have checked that the chiral and continuum extrapolations of the simulated vector meson masses are consistent with the PDG values within lattice uncertainties, which are dominated by the error of the lattice scale. where ξ ≡ M 2 π /(4πf 0 ) 2 and the exponential term is a phenomenological representation of possible finite size effects (FSEs). The results of the linear fit (4.2) are shown in figure 7 by the solid lines. In our combined fit the values of the parameters are determined by a χ 2 -minimization procedure adopting an uncorrelated χ 2 . The uncertainties on the fitting parameters do not depend on the value of the uncorrelated χ 2 , because they are obtained using the bootstrap procedure of ref. [34] (see section 3.1). This guarantees that all correlations among the lattice data points and among the fitting parameters are properly taken into account.
Averaging over the results corresponding to different fitting functions of the data either with or without the ELM procedure we get at the physical point • () input is the error coming from the uncertainties of the input parameters of the eight branches of the quark mass analysis of ref. [34]; • () disc is the uncertainty due to both discretization effects and scale setting, estimated by comparing the results obtained with and without the ELM procedure (4.1); • () FSE is the error coming from including (F s = 0) or excluding (F s = 0) the FSE correction. When FSEs are not included, all the gauge ensembles with L/a = 20 and 24 are also not included; • () chir is the error coming from including (A s 1 = 0) or excluding (A s 1 = 0) the linear term in the light-quark mass.
For each quark flavor f the e.m. correction δV (t) to the vector correlator is given by where J C µ (y) and T ν (y) are given in eqs. has been determined in ref. [29], while 1/Z m = Z P , where Z P is the RC of the pseudoscalar density evaluated in ref. [34]. For 1/Z f we use the perturbative result at leading order in α em in the MS scheme, given by [46,47] where the renormalization scale µ is taken to be equal to µ = 2 GeV, at which we consider that the renormalized quark masses in QCD and QCD+QED coincide (see ref. [29]).
Within the quenched QED approximation, which neglects the effects of the sea-quark electric charges, the correlator δV self (t) + δV exch (t) corresponds to the sum of the diagrams 8a-8b, while the correlators δV tad (t), δV PS (t) and δV S (t) represent the contributions of the diagrams 8c, 8d and 8e, respectively. In the quenched QED approximation the shift δm crit f is proportional to α em q 2 f (see for details ref. [29]). In addition one has to consider also the QED contribution to the renormalization constant of the vector current (3.23), namely where Z

(0)
A is the renormalization constant (RC) of the current in absence of QED (determined in ref. [34]) and δZ A is the O(α em ) RC. The latter can be written as where δM V and δZ V can be determined, respectively, from the "slope" and the "intercept" of the ratio δV (t)/V (t) at large time distances t min ≤ t ≤ t max (see refs. [27][28][29]). Note The labels "self", "tad + PS", "exch", "scalar" and "Z A " indicate the contributions of the diagrams 8b, 8c+8d, 8a, 8e and the one generated by the QED effect in the RC Z A of the local vector current at leading order in α em (see eq. (5.9)) with Z (fact) A = 0.9. The label "total" corresponds to the sum of all the contributions. that all the quantities δV , δZ V and δM V are proportional to α em q 2 f , which make δa had µ proportional to α 3 em q 4 f . The time dependence of the integrand function in the r.h.s. of eqs. (5.11)-(5.12) is shown in figure 9 in the case of the ETMC gauge ensemble D20.48. The contributions coming from the various diagrams of figure 8 as well as from the additional term (5.9) are determined quite precisely and are characterized by different signs. Partial cancellations among the various contributions occur in the total sum, which turns out to be smaller than each individual contributions. Thus, even a 10% uncertainty on the RC δZ A may have a larger impact on the final uncertainty of δa had µ , as it will be shown later on. The results for the strange contribution to δa had µ (<), δa had µ (>) and their sum δa had µ , obtained adopting the four choices of T data , namely: T data = (t min + 2), (t min + t max )/2, (t max − 2) and (T /2 − 4), are collected in table 4 for some of the ETMC gauge ensembles.
As in the case of the lowest-order terms a had µ (<) and a had µ (>), we find that the separation between δa had µ (<) and δa had µ (>) depends on the specific value of T data , as it should be, but their sum δa had µ is largely independent of the choice of the value of T data in the range between t min and t max within the statistical uncertainties. As in the case of the lowest-order term, the contribution δa had µ (>), which depends on the analytic representation (5.12), is significantly reduced at T data = T /2−4, where it does not exceed the statistical uncertainty of δa had µ .
In the case of the charm contribution the value of δa had µ (>) is always several orders of magnitude smaller than δa had µ (<) and the latter turns out to be the same for all the four choices of T data .    Table 4. Results for the strange contribution to δa had µ (<), δa had µ (>) and their sum δa had µ , in units of 10 −12 , obtained assuming T data = (t min + 2), (t min + t max )/2, (t max − 2) and (T /2 − 4) for the ETMC gauge ensembles A40.24, A30.32, B25.32 and D15.48. Errors are statistical only.
The precision of the lattice data can be drastically improved by forming the ratio of the e.m. correction over the lowest-order term. Therefore, in what follows we perform our analysis of the ratio δa had µ /a had µ , which is shown in figure 10. We have checked that in the case of the e.m. corrections the use of the ELM procedure (4.1) does not improve the precision of the lattice data.
It can be seen from figure 10 that the dependence on the light-quark mass m is quite mild, being driven only by sea quarks, and that the uncertainties of the data are dominated by the error on the RC δZ A , which has been taken to be the same for all the gauge ensembles used in this work (see appendix A). the structure-dependent (SD) FSEs are expected to start at order O(1/L 2 ). According to the effective field theory approach of ref. [50], one might argue that in the case of mesons with vanishing charge radius (as the ones appearing in the correlator δV (t)) the SD FSEs may start at order O(1/L 3 ). Therefore we adopt the following simple fitting function • () input is the error coming from the uncertainties of the input parameters of the eight branches of the quark mass analysis of ref. [34]; • () disc is the uncertainty due to both discretization effects and scale setting, estimated by comparing the results obtained with and without the ELM procedure (4.1); • () FSE is the error coming from including (δF s = 0) or excluding (δF s = 0) the FSE correction. When FSEs are not included, all the gauge ensembles with L/a = 24 are also not included; • () chir is the error coming from including (δA s 1 = 0) or excluding (δA s 1 = 0) the linear term in the light-quark mass.
• () Z A is the error generated by the uncertainty on the RC Z f act A (see eq. (5.8)), which turns out to be by far the dominant source of uncertainty.
Using the lowest-order results Thus, the e.m. corrections to δa s µ and δa c µ turn out to be negligible with respect to the current uncertainties of the lowest-order terms.

Conclusions
We have presented a lattice calculation of the HVP contribution of strange and charm quarks to the anomalous magnetic moment of the muon at orders O(α 2 em ) and O(α 3 em ) in the e.m. coupling.
We have employed the gauge configurations generated by the European Twisted Mass Collaboration with N f = 2 + 1 + 1 dynamical quarks at three values of the lattice spacing (a 0.062-0.089 fm) with pion masses in the range M π 210-450 MeV and with strange and charm quark masses tuned at their physical values.
In this work we have taken into account only connected diagrams, in which each quark flavor contributes separately, and a direct summation of the relevant correlators over the Euclidean time distances has been performed, adopting the local lattice version of the e.m. current operator.
As for the calculation of the e.m. corrections in the strange and charm sectors, we have adopted the RM123 approach of ref. [28], based on the expansion of the lattice path-integral in powers of the small e.m. coupling, namely α em ≈ 1%. which show that the latter ones are negligible with respect to the present uncertainties of the lowest-order terms. We stress that the current uncertainties on the e.m. corrections δa s µ and δa c µ are of the order of ∼ 60% and ∼ 40%, since they are dominated by the uncertainty on the RC Z A of the local vector current, which has been estimated through the axial Ward-Takahashi identity (WTI) derived in the presence of QED effects (see appendix A). A dedicated study aimed at the determination of the RCs of bilinear operators in presence of QED employing non-perturbative renormalization schemes, like the RI-MOM one, is expected to improve significantly the precision of the calculation of the e.m. corrections and isospin-breaking effects on a had µ . Our findings demonstrate that the expansion method of ref. [28], which has been already applied successfully to the calculation of e.m. corrections to meson masses [28,29] and to the leptonic decays of pions and kaons [30,31], works as well also in the case of the HVP contribution to the muon anomalous magnetic moment. The application of the approach presented in this work to the case of the u-and d-quark contributions is ongoing.

JHEP10(2017)157
where M (x, y) = m + i 4 r a − m crit γ 5 τ 3 δ(x, y) , with E µ (x) = e ieQAµ(x) being the QED link, A µ (x) the photon field, m the twisted bare quark mass (in QCD+QED), m crit the critical mass (in QCD+QED) and Q ≡ diag {q 1 , q 2 }. Performing the local non-singlet axial rotation where ∂ µ is the backward derivative in the µ direction and We now choose that the charges of the two quarks are the same, i.e. q 1 = q 2 = q. This implies that the isospin rotation τ + commutes with the QED link E µ (x). Consequently, the first line in eq. (A.7) vanishes, while the second and third lines can be written as a backward derivative. Thus, eq. (A.5) becomes which is conserved in the chiral limit and therefore it does not require any renormalization constant.
As is well known, the local current requires a multiplicative renormalization, given by the RC Z V [38], in order to match the 1-point split TM axial current (A.9) in the continuum limit. Thus, provided the quark charges are the same, the local version of the TM axial Ward-Takahashi identity holds as well also in the presence of electromagnetism, viz.

A.2 Determination of the RC Z V
Let's consider a pseudoscalar (PS) meson composed by the two mass-and charge-degenerate TM quarks (ψ 1 , ψ 2 ). Introducing the 2-point correlators Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.