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 $N_f = 2+1+1$ dynamical quarks at three values of the lattice spacing ($a \simeq 0.062, 0.082, 0.089$ fm) with pion masses in the range $M_\pi \simeq 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_\mu^s(\alpha_{em}^2) = (53.1 \pm 2.5) \cdot 10^{-10}$, $a_\mu^s(\alpha_{em}^3) = (-0.018 \pm 0.011) \cdot 10^{-10}$ and $a_\mu^c(\alpha_{em}^2) = (14.75 \pm 0.56) \cdot 10^{-10}$, $a_\mu^c(\alpha_{em}^3) = (-0.030 \pm 0.013) \cdot 10^{-10}$ for the strange and charm contributions, respectively.


F. Sanfilippo and S. Simula
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 one-loop diagram, a had µ (α 2 em ), as well as to the nextto-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].
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.
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 (3)(4) are found to be negligible with respect to present uncertainties.
The paper is organized as follows. In section II we introduce the basic quantities and notation.
In section III we describe the lattice calculation and give the simulation details. In section IV we present the calculation of the strange and charm contributions to the HPV at order O(α 2 em ) and in section V the corresponding e.m. corrections at order O(α 3 em ), which represent the original part of this work. Finally, section VI contains our conclusions and outlooks for future developments.

II. MASTER FORMULA
The hadronic contribution a had µ to the muon anomalous magnetic moment at order α 2 em can be related to the Euclidean space-time HVP function Π(Q 2 ) by [8][9][10] a had 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. (5) 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 cor- 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]f 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. (11) 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 with 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. (14) should be almost independent of the specific choice of the value of T data , as it will be shown later in Section IV.

A. 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 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 I [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].
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. (15)(16), 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 Fig. 1 in the case of the muon at the three values of the lattice spacing of the ETMC ensembles of Table I  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 twisted-mass 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 The variation of the lattice action with respect to a vector rotation α V (x) of the quark fields, 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 (18) 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. (23) 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. (7) and (9) and taking into account that V (t) = V (−t), one obtains In the the renormalized HVP function [Π(Q 2 ) − Π(0)] the term (Z L + Z T ) cancels out, so that Eq. (10) 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. (36) in Section III D).
C. 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 σ had (s) = 4πα 2 em s R had (s) (28) 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. (29) into Eq. (27) one obtains which exhibits a logarithmic divergence.
In the continuum the vector correlator (9) can be obtained simply by taking the Fourier transform of the spatial components of the HVP tensor (7). Choosing Q along the time direction only, one gets Using Eq. (27) one has Consequently, using the pQCD result (29) for R had (s) (including the quark mass threshold s thr = 4m 2 ) the pQCD prediction for V (t) is given by Note that e −2mt (1 + 2mt + 2m 2 t 2 ) = 1 + O(m 3 t 3 ) 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.

GeV.
A closer look to Fig. 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 (30) 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. (17) 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.

D. 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. (18) the evaluation of the vector correlator has been carried out using the following local current: where ψ f and ψ f represent two quarks with the same mass and charge, but regularized with opposite values of the Wilson r-parameter, i.e. r f = −r f . Being at maximal twist the current (36) renormalizes multiplicatively with the renormalization constant Z A of the axial current.
The choice (36) differs from the one given by Eq. (18)  and by the absence of disconnected insertions. The first point is illustrated in Fig. 3, where the contribution of the strange quark to a had µ , evaluated using either the current (36) or the connected insertion of Eq. (18), 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 (36) implies that the "purely connected" vector correlator based on the current (18) is a well defined quantity and admits the hadron decomposition necessary for having the representation (16) (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 Fig. 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   Table II.

IV. STRANGE AND CHARM CONTRIBUTIONS: LOWEST ORDER
Let's start by considering the evaluation of a had µ (<) and a had µ (>) defined in Eqs. (15)(16) for various values of the "cut" T data chosen in the range between t min and t max given in Table II. 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 III 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 (16) 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 Fig. 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 strange contribution. The errors of the data turn out to be dominated by the uncertainties of the scale setting, which are ensemble A40.24 a had µ (<) 38.03 (28) a had µ (<) 40.83 (14) 43.18 (17) The expected advantage of the use of the effective lepton mass (37) comes from the fact that the kernel function, and therefore a had µ , depends only on the lepton mass in lattice units (see Eq. (17)). Thanks to Eq. (37), 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 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 volume limits using the following simple Ansatz 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 (38) are shown in Fig. 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 III A). 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 a s,phys µ = (53.1 ± 1.6 stat+f it ± 1.5 input ± 1.3 disc ± 0.2 F SE ± 0.1 chir ) · 10 −10 , = (53.1 ± 1.6 stat+f it ± 2.0 syst ) · 10 −10 , = (53.1 ± 2.5) · 10 −10 , where • () stat+f it indicates the uncertainty induced by both the statistical errors and the fitting procedure itself; • () 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 (37); • () F SE 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.

V. STRANGE AND CHARM CONTRIBUTIONS: E.M. CORRECTIONS
Let's now turn to the e.m. corrections at leading order in α em to a had µ , which using the expansion method of Ref. [28] require the evaluation of the self-energy, exchange, tadpole, pseudoscalar and scalar insertion diagrams depicted in Fig. 8. For each quark flavor f the e.m. correction δV (t) to the vector correlator is given by with where J C µ (y) and T ν (y) are given in Eqs. (19) and (21) 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 M S scheme, given by [46,47] 1 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]).
The removal of the photon zero-mode is done according to QED L [48], i.e. the photon field A µ satisfies A µ (k 0 , k = 0) ≡ 0 for all k 0 .
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 P S (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 (36), 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 Z  [46,47], we have to add to Eq. (41) the following contribution Thus, the e.m. corrections δa had µ can be written as with (see Eqs. (15)(16)) 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 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. (51-52) is shown in Fig. 9 in the case of the ETMC gauge ensemble D20.48. The contributions coming from the various diagrams of Fig. 8 as well as from the additional term (49) 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 IV 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 (52), 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 . 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 Fig. 10. We have checked that in the case of the e.m. corrections the use of the ELM procedure (37) does not improve the precision of the lattice data.
It can be seen from Fig. 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).     is not yet available. According to the general findings of Ref. [49] the universal FSEs are expected to vanish, since they depend on the global charge of the meson states appearing in the spectral decomposition of the vector correlator δV (t). Moreover, 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 where • () stat+f it indicates the uncertainty induced by both the statistical errors and the fitting procedure itself; • () 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 (37); • () F SE 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.
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. 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%.
After extrapolation to the physical pion mass and to the continuum limit our results for a had µ are for the lowest-order contributions a s µ (α 2 em ) = (53.1 ± 2.5) · 10 −10 , 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. For an isospin doublet ψ ≡ (ψ 1 , ψ 2 ) of mass-degenerate quarks the twisted-mass (TM) action including QED is given in the physical basis at maximal twist by [36,38,39] where 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. (A7) vanishes, while the second and third lines can be written as a backward derivative. Thus, where 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 (A9) 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.
The ratio of Eq. (A21) with the lowest-order correlator C from which all the ingredients entering Eq. (A19) can be calculated.
3. Determination of the RCs Z A and Z P /Z S We now consider an isospin doublet ψ OS = ψ OS 1 , ψ OS 2 of mass-and charge-degenerate quark fields regularized using the Osterwalder-Seiler (OS) prescription [40], i.e. the same value of the Wilson r-parameter is assumed for the two quarks. At maximal twist the local axial current requires a multiplicative RC given by Z A [38]. Once renormalized the matrix elements of the TM (A10) and OS (A23) axial currents can differ only by discretization effects and therefore one has This implies where δ 0|A OS 0 |P S / 0|A OS 0 |P S (0) can be determined from the relevant axial correlators computed in the OS regularization.
Similarly, the ratio Z P /Z S of the RCs of the pseudoscalar and scalar densities can be determined by using the relation As for the e.m. corrections δZ P and δZ S one has which, however, does not allow to determine separately the two corrections δZ P and δZ S . For this reason we will not investigate Eq. (A27) numerically.

Numerical results
In order to get a first non-perturbative estimate of the RCS Z V and Z A in QCD+QED we have calculated the r.h.s. of Eqs. (A19) and (A25) using for the (bare) quark mass m the values of the strange quark masses reported in Table I for the three values of the inverse lattice coupling β of the ETMC ensembles. As in Section V, we introduce the correction factors Z (f act) V and Z (f act) A to the "naive factorization" approximation by defining where Z which cover the spread of the results given in Table V.