Electromagnetic and strong isospin-breaking corrections to the muon $g - 2$ from Lattice QCD+QED

TL;DR

This paper presents a lattice QCD calculation of electromagnetic and isospin-breaking corrections to the muon g-2, providing the most accurate estimate to date of these effects on the hadronic vacuum polarization contribution.

Contribution

The study introduces a precise lattice QCD+QED computation of isospin-breaking corrections to muon g-2, using the RM123 approach with quenched QED and multiple lattice configurations.

Findings

Calculated corrections for light, strange, and charm quarks with uncertainties.

Provided the most accurate estimate of isospin-breaking effects on muon g-2.

Results support refined theoretical predictions for muon anomalous magnetic moment.

Abstract

We present a lattice calculation of the leading-order electromagnetic and strong isospin-breaking corrections to the hadronic vacuum polarization (HVP) contribution to the anomalous magnetic moment of the muon. 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 between $\simeq 210$ and $\simeq 450$ MeV. The results are obtained adopting the RM123 approach in the quenched-QED approximation, which neglects the charges of the sea quarks. Quark disconnected diagrams are not included. After the extrapolations to the physical pion mass and to the continuum and infinite-volume limits the contributions of the light, strange and charm quarks are respectively equal to $\delta a_\mu^{\rm HVP}(ud) = 7.1 ~ (2.5) \cdot 10^{-10}$, $\delta a_\mu^{\rm HVP}(s) = -0.0053 ~ (33) \cdot 10^{-10}$ and $\delta a_\mu^{\rm HVP}(c) = 0.0182 ~ (36) \cdot 10^{-10}$. At leading order in $\alpha_{em}$ and $(m_d - m_u) / \Lambda_{QCD}$ we obtain $\delta a_\mu^{\rm HVP}(udsc) = 7.1 ~ (2.9) \cdot 10^{-10}$, which is currently the most accurate determination of the isospin-breaking corrections to $a_\mu^{\rm HVP}$.

Figures

Tables

Table 1. Table 1: Values of t m ​ i ​ n subscript 𝑡 𝑚 𝑖 𝑛 t_{min} and t m ​ a ​ x subscript 𝑡 𝑚 𝑎 𝑥 t_{max} adopted in Ref.  Giusti:2018mdh to extract the ground-state signal from the light-quark vector correlator V u ​ d ​ ( t ) superscript 𝑉 𝑢 𝑑 𝑡 V^{ud}(t) for each value of β 𝛽 \beta and of the lattice volume V / a 4 𝑉 superscript 𝑎 4 V/a^{4} for the ETMC gauge ensembles adopted in this work (see Table  4 of Appendix  A ).

$\beta$	$V/a^4$	$t_{min}/a$	$t_{max}/a$
$1.90$	${40}^3\times 80$	$12$	$22$
	${32}^3\times 64$	$12$	$22$
	${24}^3\times 48$	$12$	$20$
	${20}^3\times 48$	$12$	$20$
$1.95$	${32}^3\times 64$	$13$	$22$
	${24}^3\times 48$	$13$	$20$
$2.10$	${48}^3\times 96$	$18$	$30$

Table 2. Table 2: Values of M V u ​ d superscript subscript 𝑀 𝑉 𝑢 𝑑 M_{V}^{ud} and E n = 1 π ​ π superscript subscript 𝐸 𝑛 1 𝜋 𝜋 E_{n=1}^{\pi\pi} (see text) obtained in Ref.  Giusti:2018mdh for the three ensembles A30.32, B25.32 and D20.48, corresponding to quite similar values of the pion mass ( M π ≈ 265 ​ MeV subscript 𝑀 𝜋 265 MeV M_{\pi}\approx 265~{}\rm{MeV} ) and to different values of M π ​ L subscript 𝑀 𝜋 𝐿 M_{\pi}L (see Table  4 of Appendix  A ).

ensemble	$M_{\pi }L$	$M_V^{ud}$   (MeV)	$E_{n=1}^{\pi \pi }$   (MeV)
A30.32	3.9	843   (26)	846   (31)
B25.32	3.4	868   (25)	848   (29)
D20.48	3.0	877   (26)	839   (23)

Table 3. Table 3: Results for the contributions δ ​ a μ HVP ​ ( < ) 𝛿 superscript subscript 𝑎 𝜇 HVP \delta a_{\mu}^{\rm HVP}(<) , δ ​ a μ HVP ​ ( > ) 𝛿 superscript subscript 𝑎 𝜇 HVP \delta a_{\mu}^{\rm HVP}(>) and their sum δ ​ a μ HVP ​ ( u ​ d ) 𝛿 superscript subscript 𝑎 𝜇 HVP 𝑢 𝑑 \delta a_{\mu}^{\rm HVP}(ud) , in units of 10 − 10 superscript 10 10 10^{-10} , obtained adopting in Eqs. ( 30 - 31 ) four different choices of T d ​ a ​ t ​ a subscript 𝑇 𝑑 𝑎 𝑡 𝑎 T_{data} , namely: T d ​ a ​ t ​ a = ( t m ​ i ​ n + 2 ​ a ) subscript 𝑇 𝑑 𝑎 𝑡 𝑎 subscript 𝑡 𝑚 𝑖 𝑛 2 𝑎 T_{data}=(t_{min}+2a) , ( t m ​ i ​ n + t m ​ a ​ x ) / 2 subscript 𝑡 𝑚 𝑖 𝑛 subscript 𝑡 𝑚 𝑎 𝑥 2 (t_{min}+t_{max})/2 , ( t m ​ a ​ x − 2 ​ a ) subscript 𝑡 𝑚 𝑎 𝑥 2 𝑎 (t_{max}-2a) and ( T / 2 − 4 ​ a ) 𝑇 2 4 𝑎 (T/2-4a) for the ETMC gauge ensembles A ​ 80.24 𝐴 80.24 A80.24 , A ​ 50.32 𝐴 50.32 A50.32 , B ​ 55.32 𝐵 55.32 B55.32 and D ​ 30.48 𝐷 30.48 D30.48 (see Table  1 for the values of t m ​ i ​ n subscript 𝑡 𝑚 𝑖 𝑛 t_{min} and t m ​ a ​ x subscript 𝑡 𝑚 𝑎 𝑥 t_{max} , and Table  4 of Appendix  A for the ETMC gauge ensembles).

$T_{data}$	$\left(t_{min}+2a\right)$	$\left(t_{min}+t_{max}\right)/2$	$\left(t_{max}-2a\right)$	$\left(T/2-4a\right)$
	$0.83\left(5\right)$	$1.00\left(4\right)$	$1.13\left(5\right)$	$1.26\left(7\right)$
$\delta a_{\mu }^{\mathrm{HVP}}(>)$	$0.61\left(12\right)$	$0.43\left(10\right)$	$0.30\left(8\right)$	$0.18\left(6\right)$
$\delta a_{\mu }^{\mathrm{HVP}}\left(ud\right)$	$1.44\left(12\right)$	$1.43\left(12\right)$	$1.43\left(12\right)$	$1.44\left(12\right)$

Table 4. Table 4: Values of the simulated-quark bare masses (in lattice units), of the pion mass (in units of MeV) and of the product M π ​ L subscript 𝑀 𝜋 𝐿 M_{\pi}L for the 16 16 16 ETMC gauge ensembles with N f = 2 + 1 + 1 subscript 𝑁 𝑓 2 1 1 N_{f}=2+1+1 dynamical quarks used in this contribution (see Ref.  Carrasco:2014cwa ) and for the gauge ensemble, A40.40 added to improve the investigation of FVEs. The bare twisted masses μ σ subscript 𝜇 𝜎 \mu_{\sigma} and μ δ subscript 𝜇 𝛿 \mu_{\delta} describe the strange and charm sea doublet according to Ref.  Frezzotti:2003xj . The central values and errors of the pion mass are evaluated using the bootstrap events of the eight branches of the analysis of Ref.  Carrasco:2014cwa . The valence quarks in the pion are regularized with opposite values of the Wilson r 𝑟 r -parameter in order to guarantee that discretization effects on the pion mass are of order 𝒪 ​ ( a 2 ​ μ u ​ d ​ Λ Q ​ C ​ D ) 𝒪 superscript 𝑎 2 subscript 𝜇 𝑢 𝑑 subscript Λ 𝑄 𝐶 𝐷 {\cal{O}}(a^{2}\mu_{ud}~{}\Lambda_{QCD}) .

ensemble	$\beta$	$V/a^4$	$a{\mu }_{ud}$	$a{\mu }_{\sigma }$	$a{\mu }_{\delta }$	$N_{cfg}$	$M_{\pi }$	$M_{\pi }L$
$A40.40$	$1.90$	${40}^3\times 80$	$0.0040$	$0.15$	$0.19$	$100$	317 (12)	5.7
$A30.32$		${32}^3\times 64$	$0.0030$			$150$	275 (10)	3.9
$A40.32$			$0.0040$			$100$	316 (12)	4.5
$A50.32$			$0.0050$			$150$	350 (13)	5.0
$A40.24$		${24}^3\times 48$	$0.0040$			$150$	322 (13)	3.5
$A60.24$			$0.0060$			$150$	386 (15)	4.2
$A80.24$			$0.0080$			$150$	442 (17)	4.8
$A100.24$			$0.0100$			$150$	495 (19)	5.3
$A40.20$		${20}^3\times 48$	$0.0040$			$150$	330 (13)	3.0
$B25.32$	$1.95$	${32}^3\times 64$	$0.0025$	$0.135$	$0.170$	$150$	259 (9)	3.4
$B35.32$			$0.0035$			$150$	302 (10)	4.0
$B55.32$			$0.0055$			$150$	375 (13)	5.0
$B75.32$			$0.0075$			$80$	436 (15)	5.8
$B85.24$		${24}^3\times 48$	$0.0085$			$150$	468 (16)	4.6
$D15.48$	$2.10$	${48}^3\times 96$	$0.0015$	$0.1200$	$0.1385$	$100$	223 (6)	3.4
$D20.48$			$0.0020$			$100$	256 (7)	3.0
$D30.48$			$0.0030$			$100$	312 (8)	4.7

Table 5. Table 5: Values adopted for the coefficients Z m f ​ a ​ c ​ t superscript subscript 𝑍 𝑚 𝑓 𝑎 𝑐 𝑡 Z_{m}^{fact} [see Eq. ( 20 )] and Z A f ​ a ​ c ​ t superscript subscript 𝑍 𝐴 𝑓 𝑎 𝑐 𝑡 Z_{A}^{fact} [see Eq. ( 24 )] evaluated in Ref.  DiCarlo:2019thl for the M1 and M2 renormalization methods (see Ref.  Carrasco:2014cwa ) at the three values of β 𝛽 \beta . In Ref.  Giusti:2017jof a common value Z A f ​ a ​ c ​ t = 0.9 ​ ( 1 ) superscript subscript 𝑍 𝐴 𝑓 𝑎 𝑐 𝑡 0.9 1 Z_{A}^{fact}=0.9\,(1) , estimated through the axial Ward-Takahashi identity derived in the presence of QED effects, was adopted at all values of β 𝛽 \beta .

$\beta$	$Z_m^{fact}\left(\mathrm{M1}\right)$	$Z_A^{fact}\left(\mathrm{M1}\right)$	$Z_m^{fact}\left(\mathrm{M2}\right)$	$Z_A^{fact}\left(\mathrm{M2}\right)$
$1.90$	1.629   (41)	0.859   (15)	1.637   (14)	0.990   (9)
$1.95$	1.514   (33)	0.873   (13)	1.585   (12)	0.980   (8)
$2.10$	1.459   (17)	0.909    (6)	1.462    (6)	0.958   (3)

Equations

$$\delta a_{\mu }^{\mathrm{HVP}}(ud)$$

$$\delta a_{\mu }^{\mathrm{HVP}}(s)$$

$$\delta a_{\mu }^{\mathrm{HVP}}(c)$$

$$\delta a_{\mu }^{\mathrm{HVP}}(udsc)=7.1(2.6)(1.2{)}_{qQED+disc}\cdot 10^{-10}=7.1(2.9)\cdot 10^{-10},$$

$$a_{\mu }^{\mathrm{HVP}}=4{\alpha }_{em}^2{\int }_0^{\infty }dt{K}_{\mu }(t)V(t),$$

$$K_{\mu }(t)=\frac{4}{m_{\mu }^2}{\int }_0^{\infty }d\omega \frac{1}{\sqrt{4+{\omega }^2}}{\left(\frac{\sqrt{4+{\omega }^2}-\omega }{\sqrt{4+{\omega }^2}+\omega }\right)}^2\left[\frac{\text{mbox}cos(\omega m_{\mu }t)-1}{{\omega }^2}+\frac{1}{2}{m}_{\mu }^2{t}^2\right]$$

$$V(t)\equiv -\frac{1}{3}\sum _{i=1,2,3}\int d\vec{x}⟨J_i(\vec{x},t)J_i(0)⟩,$$

$$J_{\mu }(x)\equiv \sum _{f=u,d,s,c,...}{J}_{\mu }^f(x)=\sum _{f=u,d,s,c,...}{q}_f{\overline{\psi }}_f(x){\gamma }_{\mu }{\psi }_f(x)$$

$$a_{\mu }^{\mathrm{HVP}}(ud)=\sum _{f=u,d}[a_{\mu }^{\mathrm{HVP}}(f){]}_{(conn)}.$$

$$V(t)=V^{ud}(t)+\delta V^{ud}(t)+\dots ,$$

$$\delta V^{ud}(t)\equiv \delta V^J(t)+\delta V^T(t)+\delta V^{PS}(t)+\delta V^S(t)+\delta V^{SIB}(t)+\delta V^{{\mathcal{Z}}_A}(t)$$

$$\delta V^J(t)$$

$$\delta V^T(t)$$

$$\delta V^{PS}(t)$$

$$\delta V^S(t)$$

$$\delta V^{SIB}(t)$$

$${\mathcal{J}}_{\mu }^f(y)$$

$${\mathcal{T}}_{\nu }^f(y)$$

$$Z_m^{QCD+QED}=Z_m\left(1-4\pi {\alpha }_{em}{\mathcal{Z}}_m\right)+\mathcal{O}({\alpha }_{em}^m{\alpha }_s^n),(m>1,n\ge 0)$$

$${\mathcal{Z}}_m={\mathcal{Z}}_m^{(1)}\cdot Z_m^{fact},$$

$${\mathcal{Z}}_m^{(1)}=\frac{q_f^2}{16{\pi }^2}\left[6\text{mbox}ln(a\mu )-22.5954\right],$$

$$J_{\mu }(x)=Z_A{q}_f{\overline{\psi }}_{f^{\prime }}(x){\gamma }_{\mu }{\psi }_f(x),$$

$$Z_A^{QCD+QED}=Z_A\left(1+4\pi {\alpha }_{em}{\mathcal{Z}}_A\right)+\mathcal{O}({\alpha }_{em}^m{\alpha }_s^n),(m>1,n\ge 0)$$

$${\mathcal{Z}}_A={\mathcal{Z}}_A^{(1)}\cdot Z_A^{fact},$$

$${\mathcal{Z}}_A^{(1)}=-15.7963\frac{q_f^2}{16{\pi }^2},$$

$$\delta V^{{\mathcal{Z}}_A}(t)\equiv -0.20014\pi {\alpha }_{em}(q^4/q^2)Z_A^{fact}{V}^{ud}(t),$$

$$\delta V^{ud}(t)=\delta V^{QED}(t)+\delta V^{SIB}(t),$$

$$\delta V^{QED}(t)=\delta V^J(t)+\delta V^T(t)+\delta V^{PS}(t)+\delta V^S(t)+\delta V^{{\mathcal{Z}}_A}(t)$$

$$\delta a_{\mu }^{\mathrm{HVP}}(ud)\equiv \delta a_{\mu }^{\mathrm{HVP}}(<)+\delta a_{\mu }^{\mathrm{HVP}}(>)$$

$$\delta a_{\mu }^{\mathrm{HVP}}(<)$$

$$\delta a_{\mu }^{\mathrm{HVP}}(>)$$

$$\frac{\delta V^{ud}(t)}{V^{ud}(t)}{}_{\overrightarrow{t>>a,(T-t)>>a}}\frac{\delta Z_V^{ud}}{Z_V^{ud}}+\frac{\delta M_V^{ud}}{M_V^{ud}}{f}_{ud}(t)$$

$$f_{ud}(t)\equiv M_V^{ud}\left(\frac{T}{2}-t\right)\frac{e^{-M_V^{ud}t}-e^{-M_V^{ud}(T-t)}}{e^{-M_V^{ud}t}+e^{-M_V^{ud}(T-t)}}-1-M_V^{ud}\frac{T}{2}\approx -\left(1+M_V^{ud}t\right)$$

$$\frac{\delta a_{\mu }^{\mathrm{HVP}}(ud)}{a_{\mu }^{\mathrm{HVP}}(ud)}={\delta }_0\left[1+{\delta }_1{m}_{ud}+{\delta }_{1l}{m}_{ud}\text{mbox}ln(m_{ud})+{\delta }_2{m}_{ud}^2+D{a}^2+{\delta }_{FVE}\right],$$

Full text

I Introduction

The muon anomalous magnetic moment $a_{\mu}\equiv(g-2)/2$ is one of the most precisely determined quantities in particle physics. It is experimentally known with an accuracy of $0.54$ ppm Bennett:2006fi (BNL E821) and the current precision of the Standard Model (SM) prediction is at the level of $0.4$ ppm PDG . The discrepancy between the experimental value, $a_{\mu}^{exp}$, and the SM prediction, $a_{\mu}^{\rm SM}$, corresponds to $\simeq 3.5\div 4$ standard deviations, namely $a_{\mu}^{exp}-a_{\mu}^{\rm SM}=31.3~{}(7.7)\cdot 10^{-10}$ Jegerlehner:2017lbd , $a_{\mu}^{exp}-a_{\mu}^{\rm SM}=26.8~{}(7.6)\cdot 10^{-10}$ Davier:2017zfy and $a_{\mu}^{exp}-a_{\mu}^{\rm SM}=27.1~{}(7.3)\cdot 10^{-10}$ Keshavarzi:2018mgv .

Since the above tension may be an exciting indication of new physics (NP) beyond the SM, an intense research program is currently underway in order to achieve a significant improvement of the uncertainties. The forthcoming $g-2$ experiments at Fermilab (E989) Logashenko:2015xab and J-PARC (E34) Otani:2015lra aim at reducing the experimental uncertainty by a factor of four, down to 0.14 ppm, making the comparison of the experimental value $a_{\mu}^{exp}$ with the theoretical prediction $a_{\mu}^{\rm SM}$ one of the most stringent tests of the SM in the quest of NP effects. On the theoretical side, the main uncertainty on $a_{\mu}^{\rm SM}$ comes from hadronic contributions, related to the hadronic vacuum polarization (HVP) and light-by-light terms Jegerlehner:2017gek ; Jegerlehner:2017lbd . With the planned reduction of the experimental error, the uncertainty of the hadronic corrections will soon become the main limitation of this SM test.

The theoretical predictions for the hadronic contribution $a_{\mu}^{\rm HVP}$ 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 Davier:2010nc ; Hagiwara:2011af . An alternative approach, proposed in Refs. Lautrup:1971jf ; deRafael:1993za ; Blum:2002ii , is to compute $a_{\mu}^{\rm HVP}$ in Lattice QCD from the Euclidean correlation function of two electromagnetic (em) currents. In this respect an impressive progress in the lattice determinations of $a_{\mu}^{\rm HVP}$, which at leading order in $\alpha_{em}$ is a quantity of order ${\cal O}(\alpha^{2}_{em})$, has been achieved in the last few years Chakraborty:2014mwa ; Chakraborty:2015ugp ; Blum:2015you ; Blum:2016xpd ; Chakraborty:2016mwy ; DellaMorte:2017dyu ; Boyle:2017gzv ; Giusti:2017jof ; Chakraborty:2017tqp ; Borsanyi:2017zdw ; Blum:2018mom ; Giusti:2018mdh .

With the increasing precision of the lattice calculations, it becomes necessary to include strong and em isospin-breaking (IB) corrections, which contribute to the HVP at order ${\cal O}(\alpha^{2}_{em}(m_{d}-m_{u})/\Lambda_{QCD})$ and ${\cal O}(\alpha^{3}_{em})$, respectively. In Ref. Giusti:2017jof a lattice calculation of the IB corrections to the HVP contributions due to strange and charm quarks, $\delta a_{\mu}^{\rm HVP}(s)$ and $\delta a_{\mu}^{\rm HVP}(c)$111In the strange and charm sectors the strong IB corrections are absent at leading order in ($m_{d}-m_{u}$)., was carried out using the RM123 approach deDivitiis:2011eh ; deDivitiis:2013xla , which is based on the expansion of the path integral in powers of the mass difference ($m_{d}-m_{u}$) and of the em coupling $\alpha_{em}$. The quenched QED (qQED) approximation, which neglects the effects of sea-quark charges, was adopted and quark disconnected contractions were not included because of the large statistical fluctuations of the corresponding signals. The dominant source of uncertainty in the results of Ref. Giusti:2017jof was related to the em corrections to the renormalization constant (RC) of the local vector current, computed through the axial Ward-Takahashi identity derived in the QCD+QED theory.

In this work we present our determination of the IB corrections to the HVP contribution due to the light $u$- and $d$-quarks, $\delta a_{\mu}^{\rm HVP}(ud)$, using the same methods and lattice setup adopted in Ref. Giusti:2017jof in the case of the strange and charm contributions. A preliminary result for $\delta a_{\mu}^{\rm HVP}(ud)$ was presented in Ref. Giusti:2018vrc . Thanks to a recent nonperturbative evaluation of QCD+QED effects on the RCs of bilinear operators performed in Ref. DiCarlo:2019thl we can update the determinations of $\delta a_{\mu}^{\rm HVP}(s)$ and $\delta a_{\mu}^{\rm HVP}(c)$ made in Ref. Giusti:2017jof , obtaining a drastic improvement of the uncertainty by a factor of $\approx 3$ and $\approx 3.5$, respectively.

Within the qQED approximation and neglecting quark-disconnected diagrams the main results of the present study are:

[TABLE]

where the errors come from (statistics + fitting procedure), input parameters, chiral extrapolation, finite-volume and discretization effects. Thus, we confirm that the em corrections $\delta a_{\mu}^{\rm HVP}(s)$ and $\delta a_{\mu}^{\rm HVP}(c)$ turn out to be negligible with respect to the current uncertainties of the corresponding lowest-order terms $a_{\mu}^{\rm HVP}(s)=53.1~{}(2.5)\cdot 10^{-10}$ and $a_{\mu}^{\rm HVP}(c)=14.75~{}(0.56)\cdot 10^{-10}$ determined in Ref. Giusti:2017jof . In the case of the $u$- and $d$-quarks our finding (1) corresponds to about $1.2\%$ of the lowest-order value $a_{\mu}^{\rm HVP}(ud)=619~{}(17.8)\cdot 10^{-10}$ obtained recently in Ref. Giusti:2018mdh .

Recent calculations of the IB corrections to the HVP are: $\delta a_{\mu}^{\rm HVP}(ud)=9.0~{}(4.5)\cdot 10^{-10}$ from FNAL/HPQCD/MILC Chakraborty:2017tqp , which includes only strong IB effects, and $\delta a_{\mu}^{\rm HVP}(ud)=9.5~{}(10.2)\cdot 10^{-10}$ from RBC/UKQCD Blum:2018mom , which includes also one disconnected QED diagram. In Ref. Borsanyi:2017zdw the BMW collaboration has estimated the value $\delta a_{\mu}^{\rm HVP}(ud)=7.8~{}(5.1)\cdot 10^{-10}$ from results of the dispersive analysis of $e^{+}e^{-}$ data Jegerlehner:2017lbd . In the case of the strange contribution $\delta a_{\mu}^{\rm HVP}(s)$ RBC/UKQCD has recently obtained the result $\delta a_{\mu}^{\rm HVP}(s)=-0.0149~{}(32)\cdot 10^{-10}$ Blum:2018mom , which confirms the smallness of such contribution though it differs slightly from our finding (2).

Summing up the three contributions (1)-(3) and adding a further $\approx 15\%$ uncertainty related to the qQED approximation and to the neglect of quark-disconnected diagrams (see Section III), we get

[TABLE]

which represents the most accurate determination of the IB corrections to $a_{\mu}^{\rm HVP}$ to date.

The paper is organized as follows. In section II we describe the evaluation of the em and strong IB corrections to the light-quark HVP contribution at order ${\cal O}(\alpha^{2}_{em}(m_{d}-m_{u})/\Lambda_{QCD})$ and ${\cal O}(\alpha^{3}_{em})$ using the RM123 approach deDivitiis:2011eh ; deDivitiis:2013xla . Details of the lattice simulations are collected in the Appendix A. In section III we describe the extrapolation to the physical pion mass and to the continuum and infinite volume limits. Finally, section IV contains our conclusions and outlooks for future developments.

II Isospin-breaking corrections in the RM123 approach

We adopt the time-momentum representation for the evaluation of the HVP contribution $a^{\rm HVP}_{\mu}$ to the muon ($g-2$), namely Bernecker:2011gh

[TABLE]

where the kernel function $K_{\mu}(t)$ is given by

[TABLE]

with $m_{\mu}$ being the muon mass. In Eq. (5) the quantity $V(t)$ is the vector current-current Euclidean correlator defined as

[TABLE]

where

[TABLE]

is the em current operator with $q_{f}$ being the electric charge of the quark with flavor $f$ in units of the electron charge $e$, while $\langle...\rangle$ means the average of the $T$-product over gluon and quark fields.

We will limit ourselves to the HVP contribution of the light $u$ and $d$ quarks, indicated by $a^{\rm HVP}_{\mu}(ud)$, neglecting off-diagonal flavor terms (i.e. including quark-connected diagrams only). In this case each quark flavor $f$ contributes separately

[TABLE]

For sake of simplicity we drop the suffix $(conn)$, but it is understood that in the following we refer always to quark-connected contractions only.

In the RM123 method of Refs. deDivitiis:2011eh ; deDivitiis:2013xla the vector correlator $V(t)$ is expanded into a lowest-order contribution $V^{ud}(t)$, evaluated in isospin symmetric QCD (i.e. $m_{u}=m_{d}$ and $\alpha_{em}=0$), and a correction $\delta V^{ud}(t)$ computed at leading order in the small parameters $(m_{d}-m_{u})/\Lambda_{QCD}$ and $\alpha_{em}$:

[TABLE]

where the ellipses stand for higher order terms in $(m_{d}-m_{u})/\Lambda_{QCD}$ and $\alpha_{em}$.

The separation between the isosymmetric QCD and the IB contributions, $V^{ud}(t)$ and $\delta V^{ud}(t)$, is prescription dependent. In this work we follow Ref. Giusti:2017jof and we impose the matching condition in which the renormalized coupling and quark masses in the full theory, $\alpha_{s}$ and $m_{f}$, and in isosymmetric QCD, $\alpha_{s}^{(0)}$ and $m_{f}^{(0)}$, coincide in the $\overline{\rm MS}$ scheme at a scale of $2~{}\mbox{GeV}$. Such a prescription is known as the GRS one Gasser:2003hk .

The calculation of the IB correlator $\delta V^{ud}(t)$ requires the evaluation of the self-energy, exchange, tadpole, pseudoscalar and scalar insertion diagrams depicted in Fig. 1.

More specifically one has

[TABLE]

where $\delta V^{{\cal{Z}}_{A}}(t)$ will be described later in this Section and

[TABLE]

with ${\cal J}_{\mu}^{f}(y)$ and ${\cal T}_{\nu}^{f}(y)$ being the lattice conserved current and the tadpole operator for the quark flavor $f$, respectively222In Eqs. (17-18) the matrix $i\tau^{3}\gamma_{5}$ appears because in the twisted-mass action the Wilson term is twisted (in the so-called physical basis at maximal twist). In the case of standard Wilson fermions the matrix $i\tau^{3}\gamma_{5}$ should be replaced by the unit one.,

[TABLE]

while $\Delta_{\mu\nu}(y_{1},y_{2})$ is the photon propagator. In Eq. (14) $\delta m_{f}^{crit}$ is the em shift of the critical mass for the quark flavor $f$. In Eq. (15) the quantity ${\cal{Z}}_{m}$ is related to the em corrections to the mass RC in QCD+QED, $Z_{m}^{QCD+QED}$, as

[TABLE]

where $Z_{m}$ is the mass RC in QCD only and the product $Z_{m}\,{\cal{Z}}_{m}$ encodes the corrections at first order in $\alpha_{em}$. The quantity ${\cal{Z}}_{m}$ can be written as

[TABLE]

where ${\cal{Z}}_{m}^{(1)}$ is the pure QED contribution at leading order in $\alpha_{em}$, given in the $\overline{\rm MS}$ scheme at the renormalization scale $\mu$ by Martinelli:1982mw ; Aoki:1998ar

[TABLE]

while $Z_{m}^{fact}$ accounts for the corrections of order ${\cal{O}}(\alpha_{s}^{n})$ with $n\geq 1$ to Eq. (20). It represents the QCD correction to the “naive factorization” approximation ${\cal{Z}}_{m}={\cal{Z}}_{m}^{(1)}$ (i.e. $Z_{m}^{fact}=1$) adopted in Ref. Giusti:2017jof . Finally, Eq. (16) corresponds to the strong IB (SIB) effect (in the GRS prescription) with $m_{u}^{(0)}=m_{d}^{(0)}=m_{ud}^{(0)}$ being the renormalized light-quark mass in isosymmetric QCD.

In the numerical evaluation of the photon propagator, performed in the Feynman gauge, the photon zero-mode has been removed according to the QEDL prescription Hayakawa:2008an , i.e. the photon field $A_{\mu}$ satisfies $A_{\mu}(k_{0},\vec{k}=\vec{0})\equiv 0$ for all $k_{0}$.

In this work we make use of the same isosymmetric QCD gauge ensembles used in Ref. Giusti:2017jof , i.e. those generated by the European Twisted Mass Collaboration (ETMC) with $N_{f}=2+1+1$ dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and the charm quarks with masses close to their physical values Baron:2010bv ; Baron:2011sf . For earlier investigations of finite volume effects (FVEs) the ETMC produced three dedicated ensembles, A40.20, A40.24 and A40.32 (see Appendix A for details), which share the same light-quark mass and lattice spacing and differ only in the lattice size $L$. To improve such an investigation a further gauge ensemble, A40.40, has been generated at a larger value of the lattice size $L$.

For our maximally twisted-mass setup $\delta m_{f}^{crit}$ has been determined in Ref. Giusti:2017dmp , while $1/Z_{m}=Z_{P}$, where $Z_{P}$ is the RC of the pseudoscalar density evaluated in Ref. Carrasco:2014cwa . The coefficient $Z_{m}^{fact}$ has been recently computed in Ref. DiCarlo:2019thl in a non-perturbative framework within the RI′-MOM scheme Martinelli:1994ty .

Within the qQED approximation, which treats the dynamical quarks as electrically neutral particles, the correlator $\delta V^{J}(t)$ corresponds to the sum of the diagrams (1a)-(1b), while the correlators $\delta V^{T}(t)$ and $\delta V^{PS}(t)$ represent the contributions of the diagrams (1c) and (1d), respectively. The diagram (1e) contributes to both $\delta V^{S}(t)$ and $\delta V^{SIB}(t)$.

In our numerical simulations we have adopted the following local version of the vector current:

[TABLE]

where $\overline{\psi}_{f^{\prime}}$ and $\psi_{f}$ represent two quarks with the same mass, charge and flavor, but regularized with opposite values of the Wilson $r$-parameter (i.e. $r_{f^{\prime}}=-r_{f}$). Being at maximal twist the current (22) renormalizes multiplicatively with the RC $Z_{A}$ of the axial current. By construction the local current (22) does not generate quark-disconnected diagrams.

As discussed in Ref. Giusti:2017jof , the properties of the kernel function $K_{\mu}(t)$, given by Eq. (6), guarantee that the contact terms, generated in the HVP tensor by a local vector current, do not contribute to both $a_{\mu}^{HVP}$ and its IB correction.

Since we have adopted the renormalized vector current (22), the contribution $\delta V^{{\cal{Z}}_{A}}(t)$, appearing in Eq. (11), takes into account the em corrections to the RC $Z_{A}$ in QCD+QED, namely

[TABLE]

where $Z_{A}$ is the RC of the axial current in pure QCD (determined in Ref. Carrasco:2014cwa ), while the product $Z_{A}\,{\cal{Z}}_{A}$ encodes the corrections at first order in $\alpha_{em}$. The quantity ${\cal{Z}}_{A}$ can be written as

[TABLE]

where ${\cal Z}_{A}^{(1)}$ is the pure QED correction at leading order in $\alpha_{em}$, given by Martinelli:1982mw ; Aoki:1998ar

[TABLE]

and $Z_{A}^{fact}$ takes into account QCD corrections of order ${\cal{O}}(\alpha_{s}^{n})$ with $n\geq 1$ to Eq. (24). In this work we make use of the non-perturbative determination obtained in Ref. DiCarlo:2019thl within the RI′-MOM scheme, which improves significantly the value $Z_{A}^{fact}=0.9\,(1)$ obtained through the axial Ward-Takahashi identity in Ref. Giusti:2017jof . The values adopted for the coefficients $Z_{m}^{fact}$ and $Z_{A}^{fact}$ are collected in Table 5 of Appendix A.

Thus, the IB term $\delta V^{{\cal Z}_{A}}(t)$ is simply given by

[TABLE]

where $q^{4}/q^{2}\equiv\sum_{f=u,d}\,q_{f}^{4}/\sum_{f=u,d}\,q_{f}^{2}=17/45$ and $V^{ud}(t)$ is the lowest-order contribution of the light-quarks to the vector correlator, calculated for our lattice setup in Ref. Giusti:2018mdh .

To sum up, the IB corrections $\delta V^{ud}(t)$ can be written as the sum of two (prescription dependent) contributions as

[TABLE]

where

[TABLE]

and $\delta V^{SIB}(t)$ is given by Eq. (16).

Within the qQED approximation, where the shift $\delta m_{f}^{crit}$ is proportional to $q_{f}^{2}$ Giusti:2017dmp , and neglecting quark-disconnected diagrams the QED correlator $\delta V^{QED}(t)$ is proportional to $q^{4}\equiv\sum_{f=u,d}\,q_{f}^{4}=17/81$. Instead, the SIB correlator $\delta V^{SIB}(t)$ is proportional to $\sum_{f=u,d}\,q_{f}^{2}\,(m_{f}^{(0)}-m_{f})=(1/6)\,(m_{d}-m_{u})$. Using as inputs the experimental charged- and neutral-kaon masses the value $m_{d}-m_{u}=2.38\,(18)$ MeV was determined in Ref. Giusti:2017dmp at the physical point in the $\overline{\rm MS}(2~{}\rm GeV)$ scheme. Such a value is adopted in Eq. (16) for all gauge ensembles.

In Fig. 2 we show the dependence of both $\delta V^{QED}(t)$ and $\delta V^{SIB}(t)$ on the time distance $t$ in the case of the ETMC gauge ensemble $D20.48$ (see Appendix A).

III Results

A convenient procedure DellaMorte:2017dyu ; Giusti:2017jof ; Giusti:2018mdh consists in splitting Eq. (5) into two contributions corresponding to $0\leq t\leq T_{data}$ and $t>T_{data}$, respectively. In the first contribution the vector correlator is numerically evaluated on the lattice, while for the second contribution an analytic representation is required. If $T_{data}$ is large enough that the ground-state contribution is dominant for $t>T_{data}$ and smaller than $T/2$ in order to avoid backward signals, the IB corrections $\delta a_{\mu}^{\rm HVP}(ud)$ can be written as

[TABLE]

with

[TABLE]

where $M_{V}^{ud}$ is the ground-state mass of the lowest-order correlator $V^{ud}(t)$ and $Z_{V}^{ud}$ is the squared matrix element of the vector current between the ground-state $|V\rangle$ and the vacuum: $Z_{V}^{ud}\equiv(1/3)\sum_{i=x,y,z}$ $\sum_{f=u,d}q_{f}^{2}$ $|\langle 0|\overline{\psi}_{f}(0)\gamma_{i}\psi_{f}(0)|V\rangle|^{2}$. In Ref. Giusti:2018mdh the ground-state masses $M_{V}^{ud}$ and the matrix elements $Z_{V}^{ud}$ have been determined from a single exponential fit of $V^{ud}(t)$ using appropriate time intervals $t_{min}\leq t\leq t_{max}$, where the ground-state is dominating. For the reader’s convenience the values chosen in Ref. Giusti:2018mdh for $t_{min}$ and $t_{max}$ at each value of $\beta$ and of the lattice volume are shown in Table 1.

In Ref. Giusti:2018mdh for each ETMC gauge ensemble the lowest-order correlator $V^{ud}(t)$ was fitted at large time distances using also the two-pion finite volume spectrum. It turned out that the first two-pion energy level $E_{n=1}^{\pi\pi}$ is always close to $M_{V}^{ud}$ within the uncertainties. This is reassuring that the use of a single exponential fit in Eq. (31) reproduces properly the tail of the correlator beyond $T_{data}$. To illustrate this point we have collected in Table 2 the values of $M_{V}^{ud}$ and $E_{n=1}^{\pi\pi}$ for the three ensembles A30.32, B25.32 and D20.48, corresponding to quite similar values of the pion mass ($M_{\pi}\approx 265$ MeV) and to values of $M_{\pi}L$ ranging from $3.0$ to $3.9$ (see Table 4).

In Eq. (31) the quantities $\delta M_{V}^{ud}$ and $\delta Z_{V}^{ud}$ can be extracted respectively from the “slope” and the “intercept” of the ratio $\delta V^{ud}(t)/V^{ud}(t)$ at large time distances (see Refs. deDivitiis:2011eh ; deDivitiis:2013xla ; Giusti:2017dmp ; Giusti:2017jof ), namely

[TABLE]

where

[TABLE]

is almost a linear function of the Euclidean time $t$. This procedure is shown in Fig. 3 in the case of the gauge ensemble $D20.48$.

The time dependencies of the integrand functions $K_{\mu}(t)\,\delta V^{QED}(t)$ and $K_{\mu}(t)\,\delta V^{SIB}(t)$ are shown in Fig. 4 in the case of the ETMC gauge ensemble $B55.32$ (see Appendix A). After summation over the time distance $t$, the SIB contribution dominates over the QED one.

The results for the separate contributions $\delta a_{\mu}^{\rm HVP}(<)$ and $\delta a_{\mu}^{\rm HVP}(>)$, as well as their sum $\delta a_{\mu}^{\rm HVP}(ud)$, are obtained adopting four choices of  $T_{data}$, namely: $T_{data}=(t_{min}+2a)$, $(t_{min}+t_{max})/2$, $(t_{max}-2a)$ and $(T/2-4a)$. These results are collected in Table 3 for some of the ETMC gauge ensembles. We find that the separation between $\delta a_{\mu}^{\rm HVP}(<)$ and $\delta a_{\mu}^{\rm HVP}(>)$ depends on the specific value of $T_{data}$, as it should be, but their sum $\delta a_{\mu}^{\rm HVP}(ud)$ is independent of the specific choice of the value of $T_{data}$ within the statistical uncertainties. Note that for $T_{data}=t_{max}-2a$ the contribution $\delta a_{\mu}^{\rm HVP}(>)$, which depends on the identification of the ground-state signal, is still a significant fraction of the total value $\delta a_{\mu}^{\rm HVP}(ud)$, as it was already observed in the case of the lowest-order term $a_{\mu}^{\rm HVP}(ud)$ in Ref. Giusti:2018mdh .

All four choices of $T_{data}$ are employed in the various branches of our bootstrap analysis. The corresponding systematics is sub-dominant with respect to the other sources of uncertainties and it will not be given separately in the final error budget.

We have considered also the ratio of the IB correction $\delta a_{\mu}^{\rm HVP}(ud)$ over the leading-order term $a_{\mu}^{\rm HVP}(ud)$, which was evaluated in Ref. Giusti:2018mdh for the same gauge ensembles. The attractive feature of the ratio $\delta a_{\mu}^{\rm HVP}(ud)/a_{\mu}^{\rm HVP}(ud)$ is to be less sensitive to some of the systematic effects, in particular to the uncertainties of the scale setting. The data for $\delta a_{\mu}^{\rm HVP}(ud)$ and the ratio $\delta a_{\mu}^{\rm HVP}(ud)/a_{\mu}^{\rm HVP}(ud)$ are shown respectively in the left and right panels of Fig. 5.

It can be seen that discretization effects play a minor role, while FVEs are more relevant.

For the separate QED and SIB contributions the FVEs differ qualitatively and quantitatively. In the case of the QED data a power-law behavior in terms of the inverse lattice size $1/L$ is expected. According to the general findings of Ref. Lubicz:2016xro the universal, structure-indepedent FVEs are expected to vanish, since they depend on the global charge of the meson states appearing in the spectral decomposition of the vector correlator, while the structure-dependent (SD) FVEs start at order ${\cal{O}}(1/L^{2})$. Moreover, using the effective field theory approach of Ref. Davoudi:2014qua one may argue that in the case of mesons with vanishing charge radius (as the ones appearing in the vector correlator) the SD FVEs start at order ${\cal{O}}(1/L^{3})$ (see also Ref. Giusti:2017jof ). In the case of the SIB correlator (16), since a fixed value $m_{d}-m_{u}=2.38\,(18)$ MeV Giusti:2017dmp is adopted for all gauge ensembles, an exponential dependence in terms of the quantity $M_{\pi}L$ is expected Lin:2001ek .

In Fig. 6 the data for the QED and SIB contributions to the ratio $\delta a_{\mu}^{\rm HVP}(ud)/a_{\mu}^{\rm HVP}(ud)$ are shown in the case of the four ensembles A40.XX, which share common values of the light-quark mass and of the lattice spacing, but differ in the lattice size $L$. It can be seen that the theoretical expectations for the FVEs are consistent with the lattice data for both the QED and SIB contributions333We remind the reader that the lowest-order term $a_{\mu}^{\rm HVP}(ud)$ has nonnegligible FVEs, which are exponentially suppressed in terms of $M_{\pi}L$ Lin:2001ek (see Fig. 9 of Ref. Giusti:2018mdh )..

Since the SIB data dominate over the QED ones, the FVEs for the ratio $\delta a_{\mu}^{\rm HVP}(ud)/a_{\mu}^{\rm HVP}(ud)$ are expected to be mainly exponentially suppressed in $M_{\pi}L$.

For the combined extrapolations to the physical pion mass and to the continuum and infinite-volume limits we adopt the following fit ansatz:

[TABLE]

where the FVE term is estimated by using alternatively one of the fitting functions

[TABLE]

with $B_{0}$ and $f_{0}$ being the leading-order low-energy constants of ChPT and $\overline{M}^{2}\equiv 2B_{0}m_{ud}$. For the chiral extrapolation we consider either a quadratic ($\delta_{1l}=0$ and $\delta_{2}\neq 0$) or a logarithmic ($\delta_{1l}\neq 0$ and $\delta_{2}=0$) dependence. Half of the difference of the corresponding results extrapolated to the physical pion mass is used to estimate the systematic uncertainty due to the chiral extrapolation. Discretization effects are estimated by including $(D\neq 0)$ or excluding $(D=0)$ the term proportional to $a^{2}$ in Eq. (34). The free parameters to be determined by the fitting procedure are $\delta_{0}$, $\delta_{1}$, $\delta_{1l}$ (or $\delta_{2}$), $D$ and $F$ (or $\widehat{F}_{n}$).

In our combined fit (34) the values of the free parameters are determined by a $\chi^{2}$-minimization procedure adopting an uncorrelated $\chi^{2}$. The uncertainties on the fitting parameters do not depend on the $\chi^{2}$-value, because they are obtained by using the bootstrap samplings of Ref. Carrasco:2014cwa . This guarantees that all the correlations among the lattice data points and among the fitting parameters are properly taken into account. The quality of our fitting procedure is illustrated in Fig. 7.

At the physical pion mass and in the continuum and infinite-volume limits we get

[TABLE]

where the errors come in the order from (statistics + fitting procedure), input parameters of the eight branches of the quark mass analysis of Ref. Carrasco:2014cwa , chiral extrapolation, finite-volume and discretization effects. In Eq. (36) the uncertainty in the square brackets corresponds to the sum in quadrature of the statistical and systematic errors.

Using the leading-order result $a^{\rm HVP}_{\mu}(ud)=619.0~{}(17.8)\cdot 10^{-10}$ from Ref. Giusti:2018mdh , we obtain our determination of the leading-order IB corrections to $a^{\rm HVP}_{\mu}(ud)$, namely

[TABLE]

which comes (within the GRS prescription) from the sum of the QED contribution

[TABLE]

and of the SIB one

[TABLE]

The above results show that the IB correction (37) is dominated by the strong $SU(2)$-breaking term, which corresponds roughly to $\approx 85\%$ of $\delta a^{\rm HVP}_{\mu}(ud)$.

Our determination (37), obtained with $N_{f}=2+1+1$ dynamical flavors of sea quarks, agrees within the errors with and is more precise than both the phenomenological estimate $\delta a^{\rm HVP}_{\mu}(ud)=7.8~{}(5.1)\cdot 10^{-10}$, obtained by the BMW collaboration Borsanyi:2017zdw using results of the dispersive analysis of $e^{+}e^{-}$ data Jegerlehner:2017lbd , and the lattice determination $\delta a^{\rm HVP}_{\mu}(ud)=9.5~{}(10.2)\cdot 10^{-10}$, obtained by the RBC/UKQCD collaboration Blum:2018mom at $N_{f}=2+1$, which includes also one disconnected QED diagram. Recently, adopting $N_{f}=1+1+1+1$ simulations, the FNAL/HPQCD/MILC collaboration has found for the SIB contribution the value $\left[\delta a^{\rm HVP}_{\mu}(ud)\right]^{(SIB)}=9.0~{}(4.5)\cdot 10^{-10}$ Chakraborty:2017tqp .

Thanks to the recent nonperturbative evaluation of QCD+QED effects on the RCs of bilinear operators performed in Refs. DiCarlo:2019thl we can update the determinations of the strange $\delta a_{\mu}^{\rm HVP}(s)$ and charm $\delta a_{\mu}^{\rm HVP}(c)$ contributions to the IB effects made in Ref. Giusti:2017jof . We get

[TABLE]

to be compared with $\delta a_{\mu}^{\rm HVP}(s)=-0.018~{}(11)\cdot 10^{-10}$ and $\delta a_{\mu}^{\rm HVP}(c)=-0.030~{}(13)\cdot 10^{-10}$ given in Ref. Giusti:2017jof . The updated results confirm that the em corrections $\delta a_{\mu}^{\rm HVP}(s)$ and $\delta a_{\mu}^{\rm HVP}(c)$ are negligible with respect to the current uncertainties of the corresponding lowest-order terms $a_{\mu}^{\rm HVP}(s)=53.1~{}(2.5)\cdot 10^{-10}$ and $a_{\mu}^{\rm HVP}(c)=14.75~{}(0.56)\cdot 10^{-10}$ Giusti:2017jof . Recently Blum:2018mom in the case of the strange contribution the RBC/UKQCD collaboration has found the result $\delta a_{\mu}^{\rm HVP}(s)=-0.0149~{}(32)\cdot 10^{-10}$, which deviates from our finding (40) by $\approx 2$ standard deviations.

The sum of our three results (37), (40) and (41) yields the contribution of quark-connected diagrams to $\delta a_{\mu}^{\rm HVP}$ within the qQED approximation, namely $\delta a_{\mu}^{\rm HVP}(udsc)|_{conn}=7.1~{}(2.6)\cdot 10^{-10}$. Recently, in Ref. Blum:2018mom one QED disconnected diagram has been calculated in the case of the $u$- and $d$-quark contribution and found to be of the same order of the corresponding QED connected term. Thus, we estimate that the uncertainty related to the qQED approximation and to the neglect of quark-disconnected diagrams is approximately equal to our QED contribution (38), obtaining

[TABLE]

which represents the most accurate determination of the IB corrections to $a_{\mu}^{\rm HVP}$ to date.

Using the recent ETMC determinations of the lowest-order contributions of light, strange and charm quarks, $a^{\rm HVP}_{\mu}(ud)=619.0~{}(17.8)\cdot 10^{-10}$, $a^{\rm HVP}_{\mu}(s)=53.1~{}(2.5)\cdot 10^{-10}$ and $a^{\rm HVP}_{\mu}(c)=14.75~{}(0.56)\cdot 10^{-10}$ Giusti:2017jof ; Giusti:2018mdh , and an estimate of the lowest-order quark-disconnected diagrams, $a^{\rm HVP}_{\mu}(disc)=-12~{}(4)\cdot 10^{-10}$, obtained using the results of Refs. Borsanyi:2017zdw and Blum:2018mom , our finding (42) for the IB corrections leads to an HVP contribution to the muon ($g-2$) equal to

[TABLE]

which agrees within the errors with the recent determinations based on dispersive analyses of the experimental cross section data for $e^{+}e^{-}$ annihilation into hadrons (see Ref. Keshavarzi:2018mgv and references therein).

IV Conclusions

We have presented a lattice calculation of the isospin-breaking corrections to the HVP contribution of light quarks to the anomalous magnetic moment of the muon at order ${\cal{O}}[\alpha^{2}_{em}(m_{d}-m_{u})/\Lambda_{QCD}]$ in the light-quark mass difference and ${\cal{O}}(\alpha^{3}_{em})$ in the em coupling. We have employed the gauge configurations generated by ETMC with $N_{f}=2+1+1$ dynamical quarks at three values of the lattice spacing $(a\simeq 0.062-0.089~{}\mbox{fm})$ with pion masses in the range $M_{\pi}\simeq 210-450~{}\mbox{MeV}$ and with strange and charm quark masses tuned at their physical values determined in Ref. Carrasco:2014cwa .

The calculation of the IB corrections has been carried out adopting the RM123 approach of Refs. deDivitiis:2011eh ; deDivitiis:2013xla , which is based on the expansion of the lattice path-integral in powers of the small parameters $(m_{d}-m_{u})/\Lambda_{QCD}$ and $\alpha_{em}$, which are both of the order of ${\cal{O}}(1\%)$.

In this work we have taken into account only connected diagrams in which each quark flavor contributes separately. The leading-order em contributions to the renormalization constant of the local version of the lattice vector current, adopted in this work, have been evaluated using a recent nonperturbative calculation performed within the $RI^{\prime}$-MOM scheme in Refs. DiCarlo:2019thl . Thanks to that we have updated also the determinations of the strange $\delta a_{\mu}^{\rm HVP}(s)$ and charm $\delta a_{\mu}^{\rm HVP}(c)$ IB contributions made in Ref. Giusti:2017jof , obtaining a drastic improvement of the uncertainties.

Within the qQED approximation and neglecting quark-disconnected diagrams the main results of the present study are:

[TABLE]

Summing up the three contributions (44)-(46) and adding a further $\approx 15\%$ uncertainty related to the qQED approximation and to the neglect of quark-disconnected diagrams, we get

[TABLE]

which represents the most accurate determination of the IB corrections to $a_{\mu}^{\rm HVP}$ to date.

New QCD simulations with $N_{f}=2+1+1$ dynamical quarks close to the physical pion point Alexandrou:2018egz and the evaluation of quark-disconnected diagrams are in progress.

Acknowledgments

We gratefully acknowledge the CPU time provided by CINECA under the initiative INFN-LQCD123 on the Marconi KNL and SKL systems at CINECA (Italy). We thank S. Bacchio and B. Kostrezwa for their help in setting up the interface and the parameters for the DD$\alpha$AMG library Alexandrou:2016izb used to evaluate the quark propagators. We thank B. Kostrezwa for his help in the HMC simulations used to produce the A40.40 gauge ensemble with the tmLQCD software package Jansen:2009xp ; Abdel-Rehim:2013wba ; Deuzeman:2013xaa . G. M., V.L. and S.S. thank MIUR (Italy) for partial support under Contract No. PRIN 2015P5SBHT. G. M. thanks the partial support from ERC Ideas Advanced Grant No. 267985 “DaMeSyFla”.

Appendix A Simulation details

The ETMC gauge ensembles used in this work are the same adopted in Ref. Carrasco:2014cwa to determine the up-, down-, strange- and charm-quark masses in isospin symmetric QCD. We employ the Iwasaki action Iwasaki:1985we for gluons and the Wilson Twisted Mass Action Frezzotti:2000nk ; Frezzotti:2003ni ; Frezzotti:2003xj for sea quarks. Working at maximal twist our setup guarantees an automatic ${\cal O}(a)$-improvement Frezzotti:2003ni ; Frezzotti:2004wz .

We consider three values of the inverse bare lattice coupling $\beta$ and different lattice volumes, as shown in Table 4, where the number of configurations analyzed $(N_{cfg})$ corresponds to a separation of 20 trajectories. For earlier investigations of finite volume effects (FVEs) the ETMC had produced three dedicated ensembles, A40.20, A40.24 and A40.32, which share the same light-quark mass and lattice spacing and differ only in the lattice size $L$. To improve such an investigation a further gauge ensemble, A40.40, has been generated at a larger value of the lattice size $L$.

At each lattice spacing, different values of the sea-light-quark masses are considered. The valence- and sea-light-quark masses are always taken to be degenerate. The values of the lattice spacing in isosymmetric QCD are: $a=0.0885~{}(36),~{}0.0815~{}(30),~{}0.0619~{}(18)~{}\mbox{fm}$ at $\beta=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. Carrasco:2014cwa . 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 pseudoscalar 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 ${\cal{O}}(a^{2})$ effects, used to determine the mass RC $Z_{m}=1/Z_{P}$ in the RI′-MOM scheme.

Throughout this work the renormalized average $u/d$ quark mass $m_{ud}$ is given in the $\overline{\rm MS}$ scheme at a renormalization scale equal to $2$ GeV. We recall that, in the GRS prescription we have chosen, the renormalized average u/d quark mass in isosymmetric QCD $m_{ud}^{(0)}$ coincide with the one in QCD+QED, i.e. $m_{ud}=m_{ud}^{(0)}$, in the $\overline{\rm MS}(2\,{\rm GeV})$ scheme. At the physical pion mass ($M_{\pi}^{phys}=M_{\pi^{0}}\simeq 135\,\mbox{MeV}$) the value $m_{ud}^{phys}=3.70\,(17)\,\mbox{MeV}$ was determined in Ref. Carrasco:2014cwa , using the PDG value of the pion decay constant PDG for fixing the lattice scale.

The statistical accuracy of the meson correlator is based on the use of the so-called “one-end” stochastic method McNeile:2006bz , which includes spatial stochastic sources at a single time slice chosen randomly. In the case of the light-quark contribution we have used 160 stochastic sources (diagonal in the spin variable and dense in the color one) for each gauge configuration.

Finally, the values evaluated in Ref. DiCarlo:2019thl for the coefficients $Z_{m}^{fact}$ [see Eq. (20)] and $Z_{A}^{fact}$ [see Eq. (24)] are collected in Table 5.