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 (e.m.) 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.


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, while the current precision of the Standard Model (SM) prediction is at the level of 0.4 ppm [1]. The tension of the experimental value with the SM prediction, a exp µ − a S M µ = (28.8 ± 8.0) · 10 −10 [1], corresponds to 3.5 standard deviations and might be an exciting indication of new physics. The forthcoming g − 2 experiments at Fermilab (E989) [2] and J-PARC (E34) [3] 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 HVP correction to the one-loop diagram, a had µ (α 2 em ), as well as to the next-to-leading-order hadronic corrections, which include O(α 3 em ) contributions (see Ref. [4]).
The theoretical predictions for the hadronic contributions are traditionally obtained using dispersion relations for relating the HVP term to the experimental cross section data for e + e − annihilation into hadrons [5,6]. An alternative approach, proposed in Refs. [7][8][9], is to compute a had µ (α 2 em ) in Euclidean lattice QCD from the correlation function of two e.m. currents. In this respect an impressive progress in the lattice determinations of a had µ (α 2 em ) has been achieved in the last few years [10][11][12][13][14][15][16][17][18][19][20]. With the increasing precision of the lattice calculations, it becomes necessary to include e.m. and strong isospin breaking (IB) corrections (contributing at order O(α 3 em ) and O(α 2 em (m d − m u )), respectively) to the HVP. In this contribution we present the results of a lattice calculation of the e.m. corrections to the HVP contribution due to strange and charm quark intermediate states, obtained in Ref. [21] using the expansion method of Refs. [22,23]. Given the large statistical fluctuations, we will show only very preliminary results for the e.m. and IB corrections to the HVP contribution due to up and down quarks. For the same reason we do not have yet results for the disconnected contributions.

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 [7][8][9] where Q is the Euclidean four-momentum and f (Q 2 ) is a well-know kinematical kernel, depending also on the muon mass m µ . The HVP function [Π(Q 2 ) − Π(0)] can be determined from the vector current-current Euclidean correlator V(t) defined as where is the e.m. current with q f being the electric charge of the quark with flavor f in units of e. One gets [24] a had µ = 4α 2 wheref (t) is given byf and can be easily calculated at any value of t. In what follows we will limit ourselves to the connected contributions to a had µ . In this case each quark flavor f contributes separately. 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 0 to T /2 with T = T/a. In what follows all the overlined quantities are in lattice units. A natural procedure is to split Eq. (3) into two contributions corresponding to 0 ≤ t ≤ T data and t > T data , respectively. In the first contribution the vector correlator is directly given by the lattice data, while for the second contribution an analytic representation is required (see Refs. [16,17,19,20]). If T data is large enough that the ground-state contribution is dominant for t > T data , one can write where |V | 2 is the (squared) matrix element of the current operator between the vector ground-state and the vacuum. For each gauge ensemble the masses M V and the matrix elements Z V are extracted from a single exponential fit of the vector correlator V(t) in the range t min ≤ t ≤ t max , given explicitly in Ref. [21].

Simulation details
The ETMC gauge ensembles used in this contribution are the same adopted in Ref. [25] to determine the up, down, strange and charm quark masses. We employed the Iwasaki action for gluons and the Wilson Twisted Mass Action 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 in which the valence strange and charm quarks are regularized as Osterwalder-Seiler fermions, 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.
We considered three values of the inverse bare lattice coupling β and different lattice volumes. 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. [25]. 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.
We made use of the bootstrap samplings elaborated for the input parameters of the quark mass analysis of Ref. [25]. There, eight branches of the analysis were adopted differing in: i) 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; ii) the chiral extrapolation performed with fitting functions chosen to be either a polynomial expansion or a Chiral Perturbation Theory Ansatz in the light-quark mass; and iii) 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 .
In our numerical simulations 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 (6) renormalizes multiplicatively with the RC Z A of the axial current. The choice (6) is characterized by the absence of disconnected insertions (see Refs. [16,17,21]). The statistical accuracy of the meson correlators is based on the use of the so-called "one-end" stochastic method [26], 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 per each gauge configuration.

Lowest order
For the evaluation of the strange and charm contributions to a had µ we have adopted four choices of T data , namely: T data = (t min + 2), (t min + t max )/2, (t max − 2) and (T /2 − 4). In Ref. [21] it is shown that a had µ is almost independent of the specific choice of the value of T data . The results obtained for the strange and charm contributions to a had µ are shown by the empty markers in the lower panels of Fig. 1. We observe a mild dependence on the light-quark mass, being driven only by sea quarks, and also small residual finite size effects (FSEs) are 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 similar for all the gauge ensembles used in this contribution.
In Ref. [12] a modification of a had µ at pion masses above the physical point has been proposed in order to weaken its pion mass dependence and improve the reliability of the chiral extrapolation. Though the procedure of Ref. [12] 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 by redefining the lepton mass as The expected advantage of the use of the effective lepton mass (7) comes from the fact that the kernel function f (t), and therefore a had µ , depends only on the lepton mass in lattice units. Thanks to Eq. (7), 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 (7) in the case of the strange and charm contributions to a had µ are shown by the filled markers in Fig. 1, where the physical values for thess and cc vector masses have been taken from PDG [1] (namely, M (phys) V = 1.0195 and 3.0969 GeV, respec-tively). 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.
We have performed a combined fit for the extrapolation to the physical pion mass, the continuum and infinite volume limits using the following Ansatz where ξ ≡ M 2 π /(4π f 0 ) 2 and the exponential term is a phenomenological representation of possible FSEs. The results of the linear fit (8) are shown in Fig. 1 by the solid lines.
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 where • () stat+ f it indicates 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. [25]; • () disc is the uncertainty due to both discretization effects and scale setting, estimated by comparing the results obtained with and without the ELM procedure (7); • () FS E 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.

Electromagnetic corrections
Let's now turn to the e.m. correction δV(t) to the vector correlator at leading order in α em . Using the expansion method of Ref. [23] for each quark flavor f it can be written as where the various terms correspond to the evaluation of the self-energy, exchange, tadpole, pseudoscalar and scalar insertion diagrams (see Fig. 8 of Ref. [21]). The removal of the photon zero-mode is done according to QED L [27], i.e. the photon field A µ satisfies A µ (k 0 , k = 0) ≡ 0 for all k 0 .
In addition one has to consider the QED contribution to the RC Z A of the vector current (6): where Z (0) A is the RC in absence of QED (determined in Ref. [25]), Z (em) f /(4π) from Refs. [28], we have to add to Eq. (11) the following contribution Thus, the e.m. corrections δa had µ 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. [22,23,29]). As in the case of the lowest-order terms a had µ , we adopt for the evaluation of δa had µ the same four choices of T data . We find that δa had µ is largely independent of the choice of the value of T data within the statistical uncertainties. In the case of the e.m. corrections the use of the ELM procedure (7) does not improve the precision of the lattice data. Instead, this can be achieved by forming the ratio δa had µ /a had µ . The results for the latter are shown in Fig. 2.  It can be seen 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 f act A , which has been taken to be the same for all the gauge ensembles used in this contribution. is not yet available. According to Ref. [30] 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 correlator δV(t). Moreover, the structure-dependent (SD) FSEs are expected to start at order O(1/L 2 ). According to Ref. [31] one might argue that in the case of mesons with vanishing charge radius the SD FSEs may start at order O(1/L 3 ). Therefore we adopt the following simple fitting function δa s,c µ /a s,c µ = δA s,c 0 + δA s,c 1 m + δD s,c a 2 + δF s,c /L n (15) where the power n can be put equal to n = 2 or n = 3. In fitting our data we do not observe sensitivity to the above choices of the power n within the statistical uncertainties. At the physical pion mass and in the continuum and infinite volume limits we get δa c,phys µ a c,phys where the error budget is similar to the one described for the lowest-order results, while () Z A is the error generated by the uncertainty on the RC Z A (see Eq. (12)). The latter one is by far the dominant source of uncertainty. Using the lowest-order results (9-10) we obtain δa s,phys µ = −0.018 (11) · 10 −10 , δa c,phys µ = −0.030 (13) · 10 −10 , showing that the e.m. corrections δa s µ and δa c µ are negligible with respect to the uncertainties of the lowest-order terms. We stress that the errors appearing in Eq. (18) are dominated by the uncertainty on the RC Z A of the local vector current, estimated through the axial WTI in the presence of QED effects. A dedicated study aimed at the determination of the RCs of bilinear operators in the presence of QED employing non-perturbative renormalization schemes, like the RI-MOM one, is expected to improve the precision of the calculation of the e.m. corrections and IB effects on a had µ . Our findings demonstrate that the expansion method of Ref. [23], already applied successfully to the calculation of the e.m. corrections to meson masses [23,29] and to the leptonic decays of pions and kaons [32,33], works as well also in the case of the HVP contribution to the muon (g − 2).
In Fig. 3 we show the preliminary results for the u and d-quark contributions to a had µ and δa had µ /a had µ , based on our present limited statistics. The strong IB effect, due to the quark mass difference (m d −m u ) determined in Ref. [29], has been included in δa (u,d) µ /a (u,d) µ . It can be seen that our results for the ratio δa (u,d) µ /a (u,d) µ are in the range 0 − 1% (see also Ref. [34]). The improvement of the statistics is ongoing.