Isospin breaking corrections to meson masses and the hadronic vacuum polarization: a comparative study

TL;DR

This study computes isospin breaking and QED corrections to meson masses and the hadronic vacuum polarization using lattice QCD, comparing stochastic and perturbative methods to evaluate their precision and computational efficiency.

Contribution

It presents a comparative analysis of stochastic and perturbative approaches for including QED effects in lattice QCD calculations of meson properties.

Findings

Electromagnetic correction to muon g-2 is under 1% for up quark and 0.1% for strange quark.

Comparison of methods shows differences in precision and computational cost.

Results obtained with unphysical light quark masses, indicating scope for further physical-mass studies.

Abstract

We calculate the strong isospin breaking and QED corrections to meson masses and the hadronic vacuum polarization in an exploratory study on a $64\times24^3$ lattice with an inverse lattice spacing of $a^{-1}=1.78$ GeV and an isospin symmetric pion mass of $m_\pi=340$ MeV. We include QED in an electro-quenched setup using two different methods, a stochastic and a perturbative approach. We find that the electromagnetic correction to the leading hadronic contribution to the anomalous magnetic moment of the muon is smaller than $1\%$ for the up quark and $0.1\%$ for the strange quark, although it should be noted that this is obtained using unphysical light quark masses. In addition to the results themselves, we compare the precision which can be reached for the same computational cost using each method. Such a comparison is also made for the meson electromagnetic mass-splittings.

Tables

Table 1. Table 1: Results for the QED correction to the meson masses from the stochastic and the perturbative method. For both methods the left column shows the result in finite volume and the right panel the result in infinite volume using ( 58 ). Note, that these results have not been obtained at the physical point. The large effect on the neutral pion mass is due to a small amount of residual chiral symmetry breaking in our Domain Wall Setup (cf. discussion in section 3 and in Blum:2010ym ).

	stochastic		perturbative
	$\delta m$ $/$MeV	$\delta m^{\text{inf V}}$ $/$MeV	$\delta m$ $/$MeV	$\delta m^{\text{inf V}}$ $/$MeV
$\delta m_{{\pi }^{+}}^{\gamma }$	$3.504\pm 0.025$	$4.597\pm 0.025$	$3.459\pm 0.016$	$4.552\pm 0.016$
$\delta m_{{\pi }^0}^{\gamma }$	$1.555\pm 0.015$	$1.555\pm 0.015$	$1.538\pm 0.016$	$1.538\pm 0.016$
$\delta m_{K^{+}}^{\gamma }$	$2.722\pm 0.022$	$3.699\pm 0.022$	$2.677\pm 0.013$	$3.653\pm 0.013$
$\delta m_{K^0}^{\gamma }$	$0.547\pm 0.005$	$0.547\pm 0.005$	$0.548\pm 0.005$	$0.548\pm 0.005$

Table 2. Table 2: Comparison of pion squared mass splittings from the stochastic data and the results of a previous calculation in Blum:2010ym . Results are given in lattice units. q 1 subscript 𝑞 1 q_{1} and q 2 subscript 𝑞 2 q_{2} denote the charges of the valence quarks in units of e 𝑒 e . Both data sets use the Feynman gauge for the photon fields.

$q_1$	$q_2$	$a^2\Delta m^2$ this work	$a^2\Delta m^2$ from Blum:2010ym
2/3	2/3	$\left(5.465\pm 0.035\right)\times {10}^{-4}$	$\left(5.406\pm 0.064\right)\times {10}^{-4}$
2/3	-1/3	$\left(7.677\pm 0.052\right)\times {10}^{-4}$	$\left(7.654\pm 0.056\right)\times {10}^{-4}$
-1/3	-1/3	$\left(1.341\pm 0.009\right)\times {10}^{-4}$	$\left(1.326\pm 0.016\right)\times {10}^{-4}$

Table 3. Table 3: The HVP contribution to a μ subscript 𝑎 𝜇 a_{\mu} without QED and the QED corrections δ V ​ a μ superscript 𝛿 𝑉 subscript 𝑎 𝜇 \delta^{V}\!a_{\mu} from stochastic and perturbative data at an isospin symmetric pion mass of 340 340 340  MeV. Results have been obtained using equation ( 80 ). The last column δ V ​ a stoch − δ V ​ a pert superscript 𝛿 𝑉 superscript 𝑎 stoch superscript 𝛿 𝑉 superscript 𝑎 pert \delta^{V}\!a^{\textrm{stoch}}-\delta^{V}\!a^{\textrm{pert}} shows the correlated difference between the results from both data sets.

	$a_{\mu }^0\times {10}^{10}$	${\delta }^V{a}_{\mu }^{\text{stoch}}\times {10}^{10}$	${\delta }^V{a}_{\mu }^{\text{pert}}\times {10}^{10}$	${\delta }^V{a}_{\mu }^{\text{stoch}}-{\delta }^V{a}_{\mu }^{\text{pert}}$
$u$	$310\pm 18$	$2.6\pm 1.2$	$0.7\pm 1.2$	$1.95\pm 0.94$
$s$	$48.49\pm 0.23$	$-0.0030\pm 0.0014$	$-0.0057\pm 0.0014$	$0.0027\pm 0.0011$

Table 4. Table 4: The HVP contribution to a μ subscript 𝑎 𝜇 a_{\mu} without QED and the QED corrections from stochastic and perturbative data at an isospin symmetric pion mass of 340 340 340  MeV. Results have been obtained using Padé R 11 subscript 𝑅 11 R_{11} .

	$a_{\mu }^0\times {10}^{10}$	${\delta }^V{a}_{\mu }^{\text{stoch}}\times {10}^{10}$	${\delta }^V{a}_{\mu }^{\text{pert}}\times {10}^{10}$	${\delta }^V{a}_{\mu }^{\text{stoch}}-{\delta }^V{a}_{\mu }^{\text{pert}}$
$u$	$318\pm 11$	$0.65\pm 0.31$	$0.37\pm 0.33$	$0.27\pm 0.26$
$s$	$47.98\pm 0.25$	$-0.0030\pm 0.0012$	$-0.0049\pm 0.0011$	$0.0019\pm 0.0010$

Table 5. Table 5: Correlated difference of results obtained from using equation ( 80 ) or Padé R 11 subscript 𝑅 11 R_{11} for determining Π ^ ​ ( Q ^ 2 ) ^ Π superscript ^ 𝑄 2 \hat{\Pi}(\hat{Q}^{2})

	stoch $\times {10}^{10}$	pert $\times {10}^{10}$
$u$	$2.0\pm 1.1$	$0.3\pm 1.0$
$s$	$0.00003\pm 0.00063$	$-0.00080\pm 0.00058$

Table 6. Table 6: Results for the multiplicative renormalization of the vector current without QED Z V 0 superscript subscript 𝑍 𝑉 0 Z_{V}^{0} and the QED correction δ ​ Z V 𝛿 subscript 𝑍 𝑉 \delta Z_{V} from the perturbative and the stochastic data.

	$Z_V^0$	$\delta Z_V^{\text{stoch}}$	$\delta Z_V^{\text{pert}}$
up	$0.70209\pm 0.00083$	$-0.002674\pm 0.000043$	$-0.002756\pm 0.000044$
strange	$0.69737\pm 0.00017$	$-0.0007102\pm 0.0000016$	$-0.0007139\pm 0.0000016$

Table 7. Table 7: The QED correction to a μ subscript 𝑎 𝜇 a_{\mu} due to the QED correction to the multiplicative renormalization Z V subscript 𝑍 𝑉 Z_{V} for the local-vector current.

	${\delta }^{Z_V}{a}_{\mu }^{\text{stoch}}\times {10}^{10}$	${\delta }^{Z_V}{a}_{\mu }^{\text{pert}}\times {10}^{10}$
up	$-1.212\pm 0.052$	$-1.249\pm 0.047$
strange	$-0.04886\pm 0.00028$	$-0.04911\pm 0.00027$

Table 8. Table 8: Strong isospin breaking correction to a μ subscript 𝑎 𝜇 a_{\mu} . The middle column shows results using equation ( 80 ) to obtain Π ^ ​ ( Q ^ 2 ) ^ Π superscript ^ 𝑄 2 \hat{\Pi}(\hat{Q}^{2}) , the right column results using Padé R 11 subscript 𝑅 11 R_{11} .

	${\delta }_{\text{s}}{a}_{\mu }/q_f^2$ using (80)	${\delta }_{\text{s}}{a}_{\mu }/q_f^2$ using Padé $R_{11}$
different u and d masses	$\left(-6.1\pm 8.8\right)\times {10}^{-10}$	$\left(-6.7\pm 1.6\right)\times {10}^{-10}$
path integral expansion	$\left(-7.2\pm 7.4\right)\times {10}^{-10}$	$\left(-6.4\pm 1.7\right)\times {10}^{-10}$

Table 9. Table 9: Summary of our results for the QED correction to a μ subscript 𝑎 𝜇 a_{\mu} with an isospin symmetric pion mass of ≈ 340 absent 340 \approx 340  MeV. Results are shown for the stochastic and perturbative method.

	$a_{\mu }^0\times {10}^{10}$	${\delta }^V{a}_{\mu }\times {10}^{10}$				${\delta }^{Z_V}{a}_{\mu }\times {10}^{10}$
		stoch, (80)	stoch, $R_{11}$	pert, (80)	pert, $R_{11}$	stoch	pert
$u$	$318\left(11\right)$	$2.6\left(1.2\right)$	$0.65\left(31\right)$	$0.7\left(1.2\right)$	$0.37\left(33\right)$	$-1.212\left(52\right)$	$-1.249\left(47\right)$
$s$	$47.98\left(25\right)$	$-0.0030\left(14\right)$	$-0.0030\left(12\right)$	$-0.0057\left(14\right)$	$-0.0049\left(11\right)$	$-0.04886\left(28\right)$	$-0.04911\left(27\right)$

Table 10. Table 10: Summary of the results for the strong isospin breaking correction.

$a_{\mu }/q_f^2\times {10}^{10}$	${\delta }_s{a}_{\mu }/q_f^2\times {10}^{10}$
	diff masses, eq. (80)	diff masses, $R_{11}$	expansion, eq. (80)	expansion, $R_{11}$
$716\pm 25$	$-6.1\pm 8.8$	$-6.7\pm 1.6$	$-7.2\pm 7.4$	$-6.4\pm 1.7$

Equations

$${⟨O⟩}_0=\frac{1}{Z_0}\int \mathcal{D}[U]\mathcal{D}[\Psi ,\overline{\Psi }]O[\Psi ,\overline{\Psi },U]e^{-S_{F,0}[\Psi ,\overline{\Psi },U]}{e}^{-S_G[U]},$$

$$⟨O⟩=\frac{1}{Z}\int \mathcal{D}[U]\mathcal{D}[A]\mathcal{D}[\Psi ,\overline{\Psi }]O[\Psi ,\overline{\Psi },A,U]e^{-S_F[\Psi ,\overline{\Psi },A,U]}{e}^{-S_{\gamma }[A]}{e}^{-S_G[U]},$$

$$U_{\mu }(x)\to e^{-i{q}_{f}e{A}_{\mu }(x)}{U}_{\mu }(x),$$

$$S_{\gamma }\left[A\right]=\frac{1}{4}\sum _x\sum _{\mu ,\nu }{\left({\partial }_{\mu }{A}_{\nu }\left(x\right)-{\partial }_{\nu }{A}_{\mu }\left(x\right)\right)}^2,$$

$${\partial }_{\mu }f(x)=f\left(x+\hat{\mu }\right)-f\left(x\right).$$

$$S_{\gamma ,Feyn.}\left[A\right]=S_{\gamma }\left[A_{\mu }\right]+\frac{1}{2}\sum _x{\left(\sum _{\mu }{\partial }_{\mu }{A}_{\mu }\left(x\right)\right)}^2.$$

$$S_{\gamma ,Feyn.}\left[A\right]=-\frac{1}{2}\sum _x\sum _{\mu }{A}_{\mu }(x){\partial }^2{A}_{\mu }(x),$$

$${\partial }_{\mu }^{\ast }f(x)=f\left(x\right)-f\left(x-\hat{\mu }\right).$$

$$A_{\mu }\left(x\right)\to A_{\mu }\left(x\right)+c_{\mu },$$

$$\sum _{\vec{x}}{A}_{\mu }\left(x_0,\vec{x}\right)=0\text{for all}\mu ,x_0,$$

$$U_{\mu }(x)\to e^{-ie{q}_f{A}_{\mu }(x)}{U}_{\mu }(x).$$

$$S_{\gamma ,Feyn.}\left[A\right]=\frac{1}{2N}\sum _{k,\vec{k}\ne 0}{\hat{k}}^2\sum _{\mu }{\left|{\tilde{A}}_{\mu }\left(k\right)\right|}^2,$$

$${\hat{k}}_{\mu }=2\sin \left(\frac{k_{\mu }}{2}\right).$$

$${\left(P_C\right)}_{\mu \nu }={\delta }_{\mu \nu }-{\left|\vec{\hat{k}}\right|}^{-2}{\hat{k}}_{\mu }{\left(0,\vec{\hat{k}}\right)}_{\nu }\text{with}{\tilde{A}}_{\mu }^{\text{Coul}}(k)={\left(P_C\right)}_{\mu \nu }{\tilde{A}}_{\nu }^{\text{Feyn}}(k).$$

$$⟨O⟩={⟨O⟩}_0+\frac{1}{2}{e}^2{\frac{{\partial }^2}{\partial e^2}⟨O⟩|}_{e=0}+\mathcal{O}({\alpha }^2).$$

$$⟨O⟩={⟨O⟩}_0-\frac{(e{q}_f{)}^2}{2}{⟨O{T}_{\mu }(x)⟩}_0{\Delta }_{\mu \mu }(0)-\frac{e^2{q}_f{q}_{f^{\prime }}}{2}{⟨O{V}_{\mu }^c(x)V_{\nu }^c(y)⟩}_0{\Delta }_{\mu \nu }(x-y)+\mathcal{O}({\alpha }^2).$$

$${\Delta }_{\mu \nu }(x-y)={⟨A_{\mu }(x)A_{\nu }(y)⟩}_{\gamma }=\frac{\int \mathcal{D}[A]A_{\mu }(x)A_{\nu }(y)e^{-S_{\gamma }[A]}}{\int \mathcal{D}[A]e^{-S_{\gamma }[A]}},$$

$${\Delta }_{\mu \nu }(x-y)={\delta }_{\mu \nu }\frac{1}{N}\sum _{k,\vec{k}\ne 0}\frac{e^{ik\cdot (x-y)}}{{\hat{k}}^2},$$

$$C^{\text{exch}}(z_0)=\sum _{\mu ,\nu }\sum _{\vec{z}}\sum _{x,y}\text{Tr}\left[S^s(z,x){\Gamma }_{\nu }^c{S}^s(x,0){\gamma }_5{S}^u(0,y){\Gamma }_{\mu }^c{S}^u(y,z){\gamma }_5\right]{\Delta }_{\mu \nu }(x-y),$$

$${\Delta }_{\mu \nu }(x-y)={⟨A_{\mu }(x)A_{\nu }(y)⟩}_{\gamma }.$$

$${\Delta }_{\mu \nu }(x-y)={⟨\sum _u{\Delta }_{\mu \nu }(x-u)\eta (u){\eta }^{†}(y)⟩}_{\eta }\equiv {⟨{\tilde{\Delta }}_{\mu \nu }(x){\eta }^{†}(y)⟩}_{\eta },$$

$${⟨\eta (u){\eta }^{†}(y)⟩}_{\eta }={\delta }_{u,y}.$$

$$C^{\text{exch}}(z_0)$$

$$C^{\text{exch}}(z_0)={⟨\sum _{\mu ,\nu }\sum _{\vec{z}}\text{Tr}\left[{\Sigma }_{\mu \nu }(z,0){\gamma }_5{\Xi }_{\mu }(0,z){\gamma }_5\right]⟩}_{\eta },$$

$${\Sigma }_{\mu \nu }(z,0)$$

$${\Xi }_{\mu }(0,z)$$

$$C^{\text{self}}(z_0)$$

$$={⟨\sum _{\mu ,\nu }\sum _{\vec{z}}\text{Tr}\left[{\Lambda }_{\mu \nu }(z,0){\gamma }_5{S}^u(0,z){\gamma }_5\right]⟩}_{\eta },$$

$${\Lambda }_{\mu \nu }(z,0)=\sum _{x,y}{S}^s(z,y){\Gamma }_{\nu }^c{\tilde{\Delta }}_{\mu \nu }(y)S^s(y,x){\Gamma }_{\mu }^c{\eta }^{†}(x)S^s(x,0).$$

$$C^{\text{tad}}(z_0)$$

$$=\sum _{\mu }\sum _{\vec{z}}\text{Tr}\left[{{\rm Y}}_{\mu }(z,0){\gamma }_5{S}^u(0,z){\gamma }_5\right],$$

$${{\rm Y}}_{\mu }(z,0)=\sum _x{S}^s(z,x)T_{\mu }{\Delta }_{\mu \mu }(0)S^s(x,0).$$

Full text

1 Introduction

In recent years, lattice QCD has made remarkable progress in calculating quantities relevant to Standard Model phenomenology. Many of these calculations have reached a precision of $\lesssim 1\%$ Aoki:2016frl , e.g. the ratio $f_{K}/f_{\pi}$ of kaon and pion decay constants or the $K_{l3}$ form factor $f_{+}(0)$. Such lattice computations are usually done in the isospin symmetric limit with the masses of the up and down quarks equal ($m_{u}=m_{d}$). However, two sources of isospin breaking (IB) are present in nature. The masses of up- and down quarks are different $m_{d}\neq m_{u}$, a correction which is of the order $\mathcal{O}((m_{d}-m_{u})/\Lambda_{\textrm{QCD}})$. In addition, quarks carry an electric charge, and thus also interact electromagnetically. The latter not only applies to up and down quarks, but also to all other quark flavours. QED corrections are of $\mathcal{O}(\alpha)$, where $\alpha$ is the electromagnetic fine structure constant. Both of these isospin breaking effects are expected to be of the order of $1\%$ and thus, can no longer be neglected in applications of lattice results to phenomenology at this level of precision.

In the last few years significant progress has been made in directly including isospin breaking and QED corrections in lattice calculations. So far computations including QED on the lattice have been mainly focused on determining electromagnetic corrections to spectral quantities such as hadron masses (see e.g. Blum:2007cy ; Blum:2010ym ; Borsanyi:2013lga ; Borsanyi:2014jba ; Horsley:2015eaa ; Horsley:2015vla ; Basak:2016jnn ; Fodor:2016bgu ; deDivitiis:2013xla ; Giusti:2017dmp ). Pioneering work on the calculation of the QED correction to matrix elements has recently been published in Carrasco:2015xwa ; Lubicz:2016mpj ; Lubicz:2016xro . Another successful application of QCD$+$QED is the calculation of hadronic light-by-light scattering Blum:2014oka ; Blum:2015gfa ; Blum:2016lnc .

Two methods are commonly used to include QED in lattice QCD computations. A non-perturbative method using stochastically generated $U(1)$ gauge configurations for the photon fields was first proposed in Duncan:1996xy . We will refer to this method as the stochastic method throughout the paper (see Blum:2007cy ; Blum:2010ym ; Borsanyi:2013lga ; Borsanyi:2014jba ; Horsley:2015eaa ; Horsley:2015vla ; Basak:2016jnn ; Fodor:2016bgu for lattice QCD$+$QED calculations using the stochastic method). On the other hand, the electromagnetic coupling $\alpha$ is small in the low-energy regime, and thus, QED can be treated perturbatively. In deDivitiis:2013xla the authors proposed to expand the Euclidean path integral in orders of $\alpha$ and explicitly calculate the leading order QED corrections. We will refer to this method as the perturbative method in the following. To our knowledge, a direct comparison of results and statistical errors of both methods using the same setup and QCD gauge configurations has not yet been made.

In this paper we present an exploratory study with unphysical quark masses, in which we calculate the QED correction with both the stochastic and the perturbative methods. This allows us to directly compare results and statistical precision at the same computational cost obtained with both methods. In this study, we work in an electro-quenched setup, i.e. we consider the sea quarks as electrically neutral and mass degenerate.

Details of our strategy to include isospin breaking and QED corrections in the calculation are described in section 2. The setup of the calculation is given in section 3. In section 4 we present the results for the isospin and QED corrections to meson masses, as a starting point for comparing the perturbative and the stochastic method. In particular, we compare the statistical precision obtained from both methods and we discuss how to extract the QED correction to meson masses to consistently compare results from the stochastic and the perturbative method. In section 5 we discuss results for the QED correction to the hadronic vacuum polarization (HVP). The HVP is the leading order hadronic contribution to the anomalous magnetic moment $a_{\mu}$ of the muon. We are currently observing a $3\sigma$ Olive:2016xmw deviation between the experimentally measured value for $a_{\mu}$ and the Standard Model estimate. This has triggered increased efforts to determine the HVP in a lattice calculation (see e.g. Boyle:2011hu ; DellaMorte:2011aa ; Burger:2013jya ; Chakraborty:2014mwa ; Bali:2015msa ; Chakraborty:2016mwy ; Borsanyi:2016lpl ; DellaMorte:2017dyu ) aiming at a precision of $1\%$ to be competitive with the current most precise estimate Davier:2010nc ; Hagiwara:2011af from $e^{+}e^{-}\rightarrow$ hadrons. At this level of accuracy, isospin breaking corrections need to be included in the computation. To our knowledge, the present work constitutes the first lattice calculation of the isospin breaking corrections to the HVP. However, we wish to emphasize that this is an exploratory study at unphysical quark masses and we do not attempt to quantify finite volume effects for the QED corrections to the HVP in this study. As for the meson masses, results and statistical errors for the QED correction to the HVP calculated with the stochastic and the perturbative method are compared. Our main results and conclusions are summarized in section 6. Some preliminary results of our work have already been presented in Boyle:2016lbc .

2 Isospin Breaking on the Lattice

In the following we give details on our strategies to include isospin breaking effects. We start with a discussion of the QCD$+$QED path integral in section 2.1. The stochastic and perturbative methods to calculate the QED corrections to hadronic observables are described in sections 2.2 and 2.3, respectively. In section 2.4, we describe how we treat strong isospin breaking corrections, i.e. $m_{u}\neq m_{d}$.

2.1 Lattice QCD$+$QED Path Integral

In lattice QCD the expectation value of an observable $O$ is calculated in terms of the discretized Euclidean path integral, which is given by

[TABLE]

with quark fields $\Psi$ and $\overline{\Psi}$ and $SU(3)$ gluon fields $U$. The subscript “[math]” on $\left<O\right>_{0}$ and $S_{F,0}$ denote that these quantities are without QED.

However, quarks carry an electric charge and thus also interact electromagnetically. To account for QED effects we consider the Euclidean QED+QCD path integral

[TABLE]

with photon fields $A$. In the following, expectation values without a subscript $\left<\cdot\right>$ denote the combined QED+QCD expectation value. The observable $O$ can now, in general, also depend on the photon fields $A$ besides the quark fields $\Psi$, $\overline{\Psi}$ and the gauge fields $U$. The fermionic action $S_{F}[\Psi,\overline{\Psi},A,U]$ in (2) also contains couplings of quarks to photons and can be obtained from the action $S_{F,0}[\Psi,\overline{\Psi},U]$ by multiplying the $SU(3)$ gauge fields by appropriate $U(1)$ phases

[TABLE]

with the elementary charge $e$ and the charge $q_{f}$ of a given quark flavour, i.e. $\{q_{u},q_{d},q_{s}\}=\{2/3,-1/3,-1/3\}$. We define the non-compact photon action as

[TABLE]

with the forward derivative

[TABLE]

Here and in the following we express all quantities in units of the lattice spacing $a$. The Feynman gauge can be imposed in the photon action (4) by adding a gauge fixing term

[TABLE]

Using integration by parts, the Feynman gauge action can be written as

[TABLE]

with $\partial^{2}\equiv\sum_{\mu}\partial_{\mu}^{*}\partial_{\mu}\,$, where $\partial_{\mu}$ is the forward derivative (5) and $\partial_{\mu}^{*}$ the backward derivative defined by

[TABLE]

One important point when including QED in the lattice calculation is the treatment of the zero-mode of the photon field. This is associated with a shift symmetry of the photon action (4)

[TABLE]

which cannot be constrained by a gauge fixing condition. In our work, we choose to remove the spatial zero modes of the photon propagator on every time slice

[TABLE]

or, in momentum space $\tilde{A}_{\mu}(k_{0},\vec{k}=0)=0$. The formulation of QED resulting from this particular treatment of the zero-mode is called QEDL and was first proposed in Hayakawa:2008an . A discussion about different prescriptions of QED in a finite box with periodic boundary conditions can be found in Patella:2017fgk ; Duncan:1996xy ; Borsanyi:2014jba ; Gockeler:1989wj ; Endres:2015gda ; Lucini:2015hfa .

Throughout this paper we work in the electro-quenched approximation, i.e. when evaluating the path integral (2) we neglect QED effects in the fermion determinant $\det(D[A,U])\equiv\det(D_{0}[U])$, where $D[A,U]$ and $D_{0}[U]$ are the Dirac operators with QED and without QED, respectively. In the electro-quenched approximation effects from the electromagnetic vacuum polarization are neglected and, thus, sea quarks are electrically neutral. Effects from electro-quenching are $SU(3)$ and $1/N_{c}$ suppressed for $\mathcal{O}(\alpha)$ contributions and expected to be of the order of $\sim 10\%$ Borsanyi:2013lga of the QED correction.

The stochastic Duncan:1996xy and perturbative deDivitiis:2013xla approach that we use to include electro-quenched QED in the calculation of the path integral (2) are explained in detail in the following.

2.2 Stochastic Method

The stochastic method to include QED in lattice calculations has first been introduced in Duncan:1996xy . Since then, this method has been used in several lattice QCD $+$ QED calculations (see e.g. Blum:2007cy ; Blum:2010ym ; Borsanyi:2013lga ; Borsanyi:2014jba ; Horsley:2015eaa ; Horsley:2015vla ; Basak:2016jnn ; Fodor:2016bgu ).

In this study, we work in the electro-quenched approximation and, thus, the $U(1)$ photon gauge fields are generated independently of the $SU(3)$ gauge fields. This allows us to include QED using existing $SU(3)$ configurations. In the electro-quenched approximation sea quarks are electrically neutral. Including electromagnetic effects for the sea quarks is computationally much more expensive, since it requires either the generation of new QED$+$QCD gauge configurations, or the calculation of reweighing factors and an accompanying increase in statistical variance. Lattice calculations with dynamical QED using the stochastic method have been done in Borsanyi:2014jba ; Horsley:2015eaa ; Horsley:2015vla .

In practice one stochastically draws appropriate $U(1)$ gauge configurations for the photon fields according to the Gaussian weight $\exp(-S_{\gamma}[A])$. The new link variables are then given as the $SU(3)$ gluon gauge links multiplied by the $U(1)$ phases

[TABLE]

We define the lattice $U(1)$ photon fields at the mid-links of the lattice, i.e. we define $A_{\mu}(x)\equiv A_{\mu}(x+\hat{\mu}/2)$.

We choose to initially generate the photon fields in the Feynman gauge due to the simple structure of the action in momentum space. In momentum space the Feynman gauge action (7) is given by

[TABLE]

where $\tilde{A}_{\mu}(k)$ is the photon field in momentum space, $N$ is the total number of lattice points and the lattice momentum $\hat{k}_{\mu}$ is given by

[TABLE]

The sum over $k,\vec{k}\neq 0$ in equation (12) indicates the removal of all spatial zero modes.

Equation (12) implies, that all components $\tilde{A}_{\mu}(k)$ of the photon field can be drawn independently of each other from a Gaussian distribution with variance $2N/\hat{k}^{2}$.

To check for gauge invariance in our calculation, we use photon fields in the Feynman and the Coulomb gauge. A Feynman gauge photon field can be transformed into the Coulomb gauge by using an appropriate projector Borsanyi:2014jba

[TABLE]

After generating the photon field in momentum space, it is converted to position space using a Fast Fourier Transform (FFT)111http://www.fftw.org/. Once the photon field configurations are transformed into position space, they are multiplied with the $SU(3)$ gauge links according to equation (11). The calculation of hadronic observables then proceeds as in the case without QED, but using the combined QED$+$QCD gauge configurations. With the stochastic method QED corrections are calculated to all orders in $\alpha$ at once albeit in the electro-quenched approximation.

Although the leading order QED corrections are of $\mathcal{O}(e^{2})$, the statistical noise contains contributions at $\mathcal{O}(e)$, which would vanish in the limit of infinitely many QED configurations because of charge conjugation invariance. However, this $\mathcal{O}(e)$ noise can be exactly removed on every gauge configuration by averaging over calculations using $+e$ and $-e$ Blum:2007cy . Since we are interested in QED corrections to hadronic quantities, we calculate correlation functions once without QED ($e=0$) and once with QED, while averaging over $+e$ and $-e$. Thus, the stochastic method requires $3$ inversions per quark flavour and source position ($e=0$, $+e$ and $-e$).

2.3 Perturbative Method

In addition to the stochastic method to include QED in our lattice calculation, we use a perturbative method, adopting the approach developed in deDivitiis:2013xla . We will summarize this method in section 2.3.1 and give details on our strategy to calculate the required correlation functions in section 2.3.2. The perturbative method has been used in Carrasco:2015xwa ; Lubicz:2016mpj to determine the QED corrections to matrix elements.

2.3.1 Introduction

Since the electromagnetic coupling $\alpha$ is small in the low-energy regime, QED can be treated perturbatively. This is done by expanding the path integral (2) as a series in the electromagnetic coupling

[TABLE]

At leading order, $\mathcal{O}(\alpha)$, one finds contributions with either two insertions of the conserved vector current $V^{c}_{\mu}$ or one insertion of the tadpole operator $T_{\mu}$ deDivitiis:2013xla

[TABLE]

Note, that equation (16) is only valid, when the operator $O$ does not depend on the electromagnetic coupling $e$. If the operator itself depends on $e$, this has to be taken into account, when expanding the path integral222This is not relevant for the QED correction to meson masses. However, for the HVP we use a setup with a conserved vector current at the sink, which in the lattice discretized theory contains a link variable, and thus, including QED, depends on $e$. A more detailed discussion can be found in section 5.2..

The conserved vector current and the tadpole operator for the Domain Wall fermion action used in this work are given in (105) and (106), respectively. The photon propagator $\Delta_{\mu\nu}(x-y)$ in equation (16) is given as

[TABLE]

with $\mu,\nu=1,\ldots,4$. In the Feynman gauge the photon propagator is given by

[TABLE]

where we subtract all spatial zero modes, i.e. we use the QEDL formulation Hayakawa:2008an as in the stochastic approach above. In addition, to numerically check for gauge invariance of the observables studied in this work, we use the Coulomb gauge. The photon propagator in the Coulomb gauge is given in the appendix in equation (127).

For mesonic two-point functions one obtains from equation (16) at leading order in $\alpha$ three different types of quark-connected Wick contractions: a photon exchange diagram, a quark self-energy diagram and a tadpole diagram. These diagrams are shown in figure 1.

We do not include any quark-disconnected diagrams in our study. In particular, we neglect diagrams that correspond to photons coupling to sea quarks. These diagrams would originate from an expansion of the fermion determinant, however, we work in the electro-quenched approximation where QED effects for the fermion determinant are neglected. In the perturbative method working in unquenched QED is possible by additionally calculating the appropriate quark-disconnected diagrams. For the stochastic method unquenched QED requires the generation of combined QCD+QED gauge configurations at substantial extra cost.

2.3.2 Numerical Calculation

To illustrate how we calculate the diagrams shown in figure 1, we now consider as an example the photon exchange diagram for a charged kaon. The corresponding correlation function is given by

[TABLE]

where $\Gamma^{c}_{\mu}$ denotes a conserved vector current insertion $V^{c}_{\mu}(x)\equiv\overline{\Psi}(x)\Gamma^{c}_{\mu}\Psi(x)$ and $S^{f}(0,x)$ is the propagator from [math] to $x$ for a quark of flavour $f$.

We calculate the correlation functions such as (19) using sequential propagators. For this, the photon propagator has to be factorized into a factor that depends only on the position $x$ of one of the photon vertices and another factor that depends only on the position $y$ of the other photon vertex, such that sequential sources with insertions of the conserved vector current and a respective factor of $\Delta_{\mu\nu}(x-y)$ at $x$ or at $y$ can be constructed. This factorization can be achieved by inserting sets of stochastic sources in the photon propagator. In this work, we will do this in two different ways, which lead to different numerical costs and different statistical errors. One possibility is to use the same stochastic source for all Lorentz indices $\mu,\nu$ of the photon propagator and to calculate sequential propagators for every combination of ${\mu,\nu}$ separately. We will call this method single-$\mu$ insertion in the following. On the other hand, one can use four different stochastic sources deDivitiis:2013xla – one for every Lorentz index $\mu={1,2,3,4}$ – and to include the sum over $\mu$ or $\nu$ in equation (19) already in the sequential source. We will call this method summed-$\mu$ insertion in the following. Both methods will be illustrated below.

We note, that it is also possible to use the stochastic photon fields generated for the stochastic method as an insertion at $x$ and $y$ Giusti:2017dmp by using

[TABLE]

However, this is simply the exact $\mathcal{O}(\alpha)$-truncation of the stochastic method. Higher order $\mathcal{O}(\alpha^{2})$ effects that differentiate the the two methods are, as we argue later, small and barely significant at this level of precision. Thus, we consider the setup of Giusti:2017dmp effectively identical to the stochastic approach and we did not perform a dedicated calculation to reproduce it.

Single-$\mu$ Insertion

For the numerical calculation of the correlation functions such as (19), we rewrite the photon propagator as

[TABLE]

with stochastic sources $\eta$ that fulfil the condition

[TABLE]

Here, we choose complex $\mathbb{Z}_{2}$ noise sources $\eta(x)$, that have randomly picked entries from $\left\{\frac{1}{2}\left(\pm 1\pm i\right)\right\}$ for every lattice site. The insertion of the set of stochastic sources in (21) allows to factorize the photon propagator with a factor $\tilde{\Delta}_{\mu\nu}(x)$ that only depends on the position $x$ of one of the photon vertices and another factor $\eta^{\dagger}(y)$ that only depends on the position $y$ of the other photon vertex. We calculate $\tilde{\Delta}_{\mu\nu}(x)$ using a Fast Fourier Transform. The correlation function (19) for the photon exchange for a charged kaon can now be written as

[TABLE]

We construct this correlation function $C(z_{0})$ using sequential propagators with insertions of the conserved vector current and either $\tilde{\Delta}_{\mu\nu}(x)$ or $\eta^{\dagger}(y)$

[TABLE]

with the sequential propagators

[TABLE]

To build correlation functions of the type quark self-energy, we use appropriate double sequential propagators, e.g. the quark self-energy diagram for a charged kaon with the photon attached to the $s$ quark is calculated as

[TABLE]

with the sequential propagator

[TABLE]

The tadpole diagrams can be constructed from a sequential propagator with an insertion of the tadpole operator multiplied with the tadpole value $\Delta_{\mu\mu}(0)$ of the photon propagator, e.g. the tadpole diagram for a charged kaon with the photon attached to the $s$ quark is calculated as

[TABLE]

with the sequential propagator

[TABLE]

The tadpole value $\Delta_{\mu\mu}(0)$ of the photon propagator can be calculated exactly for a given lattice size. All the required building blocks to construct the quark-connected diagrams for the $\mathcal{O}(\alpha)$ QED correction diagrams for mesonic two-point functions are shown in figure 2.

The evaluation of the correlation functions shown in figure 1 using sequential propagators as described above requires a total of $17$ inversions per quark flavour and source position if the photon propagator is in the Feynman gauge, where only diagonal terms $\mu=\nu$ contribute (cf. equation (18)). These $17$ inversions are split as follows: $1$ inversion for the point-to-all propagator, $4$ sequential inversions with an insertion of $\Gamma^{c}_{\nu}\,\tilde{\Delta}_{\mu\nu}(y)$ (for the Feynman gauge only $\mu=\nu$ is required), $4$ sequential inversions with insertion of $\Gamma^{c}_{\mu}\,\eta^{\dagger}(x)$, $4$ additional inversions to obtain the double sequential propagators and $4$ sequential inversions for the tadpole using $T_{\mu}\Delta_{\mu\mu}(0)$ as insertion. If one uses a different gauge (e.g. Coulomb gauge) where also off-diagonal terms $\mu\neq\nu$ contribute, more inversions are required. Thus, in terms of numerical cost, the Feynman gauge is favourable for the perturbative approach with this setup to calculate sequential propagators. However, we also calculate the $\mathcal{O}(\alpha)$ QED correction using the Coulomb gauge on a subset of the statistics to check for gauge invariance.

We note, that the insertion of the photon propagator can be done using stochastic sources at both vertices by

[TABLE]

with two sets of stochastic sources $\eta$ and $\zeta$. The correlation functions that determine the QED corrections to the mesonic two-point functions are then calculated using sequential sources with appropriate insertions of either $\eta^{\dagger}$ or $\zeta^{\dagger}$. The photon propagator $\Delta_{\mu\nu}$ is included as

[TABLE]

which, for a given combination of stochastic sources $\eta^{\dagger}$, $\zeta^{\dagger}$ and $\mu$, $\nu$, is simply an overall numerical factor which multiplies the remainder of the correlation function after all the quark contractions have been calculated. This allows us to study e.g. different gauges or QED prescriptions without having to calculate new quark contractions, and thus, new quark inversions. However, for the setup that we use in this exploratory study, this resulted in a significantly worse noise-to-signal ratio for the QED corrections. Inserting the photon propagator stochastically at both vertices we found the statistical error to be $\approx 30$ times larger for the photon exchange diagram and $\approx 60$ times larger for the quark self-energy diagram compared to using only one set of stochastic sources. Thus, for the study presented here, we decided to use only one set of stochastic sources at one of the photon vertices.

Summed-$\mu$ Insertion

The number of inversions required for the construction of the diagrams shown in figure 1 can be substantially reduced by using different stochastic sources for the $4$ Lorentz indices of the photon propagator deDivitiis:2013xla ; Giusti:2017dmp . We start by rewriting the photon propagator as

[TABLE]

with stochastic sources

[TABLE]

The photon exchange diagram (19) can now be written as

[TABLE]

with the sequential propagators

[TABLE]

Each of these sequential propagators can be calculated with a single inversion using a sequential source with an insertion of either $\sum_{\nu}\Gamma^{c}_{\nu}\,\hat{\Delta}_{\nu}(x)$ or $\sum_{\mu}\Gamma^{c}_{\mu}\,\xi^{\dagger}_{\mu}$. The sequential propagators (39) and (40) are depicted in figure 3.

The calculation of the diagrams in figure 1 requires a total of $5$ inversions per quark flavour and source position with the summed-$\mu$ insertion, $1$ for the point-to-all propagator, $1+1$ for the two sequential propagators (39) and (40), $1$ additional inversion for the double sequential propagator for the quark self-energy diagram and $1$ inversion for the tadpole diagram using $\sum_{\mu}T_{\mu}\Delta_{\mu\mu}(0)$ as sequential insertion. Thus, the summed-$\mu$ insertion method is cheaper in computational cost compared to the single-$\mu$ insertion. However, the statistical error is expected to be larger, since the unwanted combinations, e.g. $\mu\neq\nu$ for the Feynman gauge, will contribute to the statistical noise.

In this study, we use both, the single- and summed-$\mu$ insertion methods, and compare the statistical precision with the stochastic method to include QED in the lattice calculation.

2.4 Strong Isospin Breaking

Even in the absence of QED, i.e. in pure QCD, isospin symmetry is broken by the different bare masses of up and down quarks. In this work, we use two different strategies to account for effects from the strong isospin breaking, one by putting different values for the valence up- and down-quark masses, and one by expanding the Euclidean path integral in the quark mass deDivitiis:2011eh .

Strong isospin breaking can be treated in a lattice calculation by simply using different values for up- and down-quark masses. In Fodor:2016bgu the up- and down-quark mass difference has been determined in the $\overline{\textrm{MS}}$ scheme at $2$ GeV as

[TABLE]

In section 3 we specify the values we choose for the bare up- and down-quark mass for the setup used in this work to approximately reproduce the physical quark mass difference (41). We include strong isospin breaking in a quenched setup, i.e. keeping the isospin symmetric sea-quark masses, to avoid having to generate new gauge configurations.

In addition, we use a strategy proposed in deDivitiis:2011eh to account for strong isospin corrections. The idea is, to expand the path integral around the isospin symmetric light quark mass $\hat{m}$

[TABLE]

where $m_{f}$ is either the mass of the down quark ($f=d$) or the up quark ($f=u$). In this way, one explicitly calculates the leading isospin breaking correction $\mathcal{O}(m_{f}-\hat{m})$. Evaluating the derivative in equation (42) one finds

[TABLE]

with the scalar current

[TABLE]

for a quark field $\psi_{f}$ of flavour $f$. A detailed derivation of equation (43) for the Domain Wall Fermion action used in this study can be found in appendix B.3.

For a mesonic two-point function one finds at $\mathcal{O}(m_{f}-\hat{m})$ one type of quark-connected contribution, which is shown in figure 4. Note, that we do not include strong isospin breaking effects for the sea quarks.

To illustrate how we calculate the diagram shown in figure 4, we consider as an example the strong isospin correction for a kaon , which is, according to equation (43), determined by a correlation function of the form

[TABLE]

with a light quark propagator $S^{l}$ using the isospin symmetric quark mass $\hat{m}$. We calculate (45) using a sequential propagator

[TABLE]

Both approaches for the inclusion of strong isospin breaking effects used in this study, are equal in terms of computational cost, i.e. number of inversions, since it requires either one inversion per gauge configuration and source position using a different quark mass for the down quark, or one inversion per gauge configuration and source position to calculate the sequential propagator $\Omega(z,0)$ using the isospin symmetric quark mass. However, using the expansion of the path integral in the mass is more flexible, since the deviation of the quark mass from the isospin symmetric mass $(m_{f}-\hat{m})$ is a free parameter, which is multiplied to the correlation function after all quark contractions have been computed. This allows for tuning the quark masses a posteriori, e.g. for fixing hadron masses to their physical value.

3 Computational Setup

For this study we use a $64\times 24^{3}$ lattice with $N_{f}=2+1$ dynamical flavours of Shamir Domain Wall Fermions Kaplan:1992bt ; Shamir:1993zy with a Domain Wall height of $M_{5}=1.8$ and $L_{s}=16$, where $L_{s}$ is the length of the fifth dimension. For further details see Allton:2008pn ; Aoki:2010dy . This gauge ensemble has been generated by the RBC/UKQCD collaboration using the Iwasaki gauge action Iwasaki:1984cj ; Iwasaki:1985we . The inverse lattice spacing of this ensemble has been determined without QED as $a^{-1}=1.78$ GeV Boyle:2017jwu . The bare sea quark masses are $am_{l}=0.005$ and $am_{s}=0.04$ for light and strange quarks, respectively. With such quark masses the isospin symmetric pion mass on this QCD gauge ensemble is $m_{\pi}=340$ MeV, thus, in this work we do not calculate at physical quark masses, even in the absence of QED.

In the present study, we use different values for the valence up- and down-quark masses. While keeping the valence up-quark at the same mass as the light quarks in the sea, we choose the valence down quark mass as $am_{d}=0.005915$. Using the results from Blum:2014tka to convert the bare mass difference $a(m_{d}-m_{u})=0.000915$ to $\overline{\textrm{MS}}$, we find $m^{R}_{d}-m^{R}_{u}=2.4$ MeV for $\overline{\textrm{MS}}$ at $2$ GeV, and thus we reproduce the physical mass difference given in Fodor:2016bgu (cf. equation (41)).

This choice ignores any QED effects in the renormalization of the quark mass. While this is acceptable for the comparative study presented here, more work is needed when aiming at physical predictions, see e.g. Borsanyi:2014jba ; Horsley:2015vla . For instance, we know that there is a small additive correction to the quark mass under renormalization in our setup. The correction can be quantified in terms of the residual mass which Blum:2010ym determined to be $m_{\rm res}\approx 0.003$ on the above ensemble (see Blum:2010ym for details). The residual mass is defined such that in the chiral limit $m_{f}=-m_{\rm res}$ and it receives additional contributions in QCD+QED which are of order $\mathcal{O}(\alpha m_{\rm res})$. Moreover, the multiplicative renormalization of the quark mass will receive QED contributions at $\mathcal{O}(\alpha)$ which have not been taken into account here.

A consequence of these simplifications in our choice of parameters is that the neutral pion splitting in the chiral limit does not vanish for finite $L_{s}$ Blum:2010ym and indeed, in this work we find the neutral pion mass shift due to QED to be sizeable. Since we are only interested in a comparative study of approaches to Lattice QCD+QED no attempt has been made to correct for this effect. Note that this effect is much more severe for lattice quark actions not obeying chiral symmetry such as Wilson fermions Portelli:2010yn .

For the bare valence strange-quark mass we use $am_{s}=0.03224$ Blum:2014tka , which, without QED, corresponds to the physical strange quark mass.

Working with physical quark masses, requires to tune the quark masses to their physical values including QED. This could be done, for example, by tuning the up-, down- and strange-quark masses until the masses of charged pion and neutral and charged kaons agree with their experimentally measured values. In addition, this requires the determination of the lattice spacing including QED, which could be done by fixing another hadron mass to its physical value, e.g. the $\Omega$-baryon. However, since this is an exploratory study and mainly focused on the comparison of the stochastic and perturbative method for including QED, we have not retuned any of the quark masses in the presence of QED.

We use $87$ QCD gauge configurations and $16$ source positions with $\mathbb{Z}_{2}$ wall sources Foster:1998vw ; McNeile:2006bz ; Boyle:2008rh for the quark propagators. For the stochastic method we use one $U(1)$ QED configuration per QCD gauge configuration. For the perturbative method we use one $\mathbb{Z}_{2}$ noise for the insertion of the photon propagator per QCD gauge configuration and source position for the single-$\mu$ insertion method and one $\mathbb{Z}_{2}$ noise for every Lorentz index for the summed-$\mu$ insertion.

4 Isospin Breaking Corrections to Meson Masses

As a starting point for comparing results from the stochastic and the perturbative method we calculate the isospin breaking corrections to meson masses. Several other calculations of QED corrections to meson masses exist even at, or extrapolated to, the physical point, see e.g. Duncan:1996xy ; Blum:2007cy ; Blum:2010ym ; Borsanyi:2013lga ; Borsanyi:2014jba ; deDivitiis:2013xla ; Horsley:2015eaa ; Horsley:2015vla ; Basak:2016jnn ; Fodor:2016bgu ; Giusti:2017dmp . In this work, we use an exploratory setup with one gauge ensemble at non-physical quark masses. However, for the first time, we directly compare results from the stochastic and perturbative methods. We also explain that the QED correction to meson masses has to be determined in different ways for the stochastic and the perturbative data, to obtain results, which can properly be compared to each other.

4.1 Extraction of the QED Correction to the Effective Mass

The two-point correlation function for a pseudoscalar meson interpolation operator $\left(\overline{\psi}_{f}\,\gamma_{5}\,\psi_{f^{\prime}}\right)$ with quark flavours $f$ and $f^{\prime}$ and vanishing spatial momentum $\vec{p}=0$, which is created at [math] and annihilated at $x$, is given by

[TABLE]

Such a two-point correlation function has the following time-dependence for large Euclidean times, where excited-state contributions are suppressed

[TABLE]

for a lattice with time extend $T$ and periodic boundary conditions. The parameter $m$ that determines the leading exponential decay of (48) is the mass of the ground state meson, whereas excited-state contributions are exponentially suppressed. In this study we are interested in the mass $m$ of the ground state. A common method to determine the mass of the ground state meson from a two-point correlation function $C(t)$ is to calculate an effective mass. In this work we use the definition of the effective mass, where one solves

[TABLE]

for $m_{\mathrm{eff}}$ at every $t$.

In the following, we discuss how to determine the QED correction to the effective mass. The effective mass including QED is given by the effective mass $m^{0}_{\textrm{eff}}$ without QED plus the QED correction $\delta m_{\textrm{eff}}$

[TABLE]

For the data from the stochastic approach the two-point function including QED contains corrections to all orders in $\alpha$ and has the form (48) with $A=A_{0}+\delta A$ and $m=m_{0}+\delta m$. Thus, the QED correction to the effective mass can be obtained by determining the effective mass according to equation (49) once for the two-point function with QED and once for the two-point function without QED and taking their difference

[TABLE]

In the following we refer to this method to extract the QED correction to the effective mass as the cosh-mass method, which is the appropriate method to extract the QED correction using the stochastic data.

On the other hand, the two-point function including QED can be expanded deDivitiis:2013xla (neglecting the backwards propagating meson for simplicity)

[TABLE]

In the perturbative approach, one explicitly only calculates the QED correction to the two-point function which are of $\mathcal{O}(\alpha)$. Keeping only terms which are of order $\alpha$ in equation (52) one finds

[TABLE]

for the QED correction $\delta C(t)$ from the perturbative data. Equation (53) implies, that the QED correction to the effective mass can be defined from the ratio of the QED correction $\delta C(t)$ to the two-point function and the two-point function $C_{0}(t)$ without QED

[TABLE]

Equation (54) can be extended to include the effects of the periodic boundary conditions

[TABLE]

using $m_{0}$ from a determination from the two-point function without QED as an input. In the following we refer to this method to extract the QED correction to the effective mass as the ratio method, which is the appropriate method to extract the QED correction using the perturbative data.

When using the ratio method for the stochastic data, one has to take into account, that the QED correction to the two-point function includes QED corrections to all orders in $\alpha$. Keeping also higher order terms in the expansion (52) one finds

[TABLE]

for the ratio method from the stochastic data. The underlined terms in (56) are included in the stochastic data, but not in the perturbative data. Thus, one expects the QED correction to the effective mass extracted using the cosh-mass method (51) and the ratio method (55) from the stochastic data to differ by

[TABLE]

We indeed find this difference in our data as illustrated in figure 5. The plot on the left-hand side shows the QED correction to the effective mass from the stochastic data extracted with the cosh-mass method (blue squares) and the ratio method (purple circles). We find a significant difference between the results from both extraction methods. The correlated difference is plotted on the right-hand side of figure 5. We can numerically confirm, that this difference is given by equation (57) as expected.

Thus, in the following we determine the QED correction to meson masses using the cosh-mass method for the stochastic data and the ratio method for the perturbative data.

4.2 QED Correction to Meson Masses

In the following we show results for the QED correction to meson masses. In subsection 4.2.1 the QED corrections to meson masses are determined and results from the perturbative and the stochastic method are compared. In subsection 4.2.2 we compare the statistical errors on the results from both methods.

4.2.1 Results

The left-hand side of figure 6 shows the QED correction to the effective mass of a charged kaon using the Feynman gauge for the photon fields. The red squares show results from the perturbative data using the ratio method (55) to extract the QED correction and the blue circles are results from the stochastic data using the cosh-mass method (51). For the perturbative data the results shown have been calculated using the single-$\mu$ insertion, which, for the same amount of statistics, gives a smaller statistical error than the summed-$\mu$ insertion (see section 4.2.2 for a detailed comparision of statistical errors).

The plot on the right-hand side of figure 6 shows the correlated difference between the stochastic and perturbative data. Both datasets are correlated since they have been calculated on the same QCD gauge configurations and the same source positions with the same $\mathbb{Z}_{2}$ wall sources for the quark propagators. Statistical errors are estimated using the bootstrap resampling method.

We find the correlated difference between both datasets to be non-zero at the level of $\approx 1.5\sigma$ and of the order of $1\%$ of the QED correction itself, which can be attributed to $\mathcal{O}(\alpha^{2})$ effects. To check this, we have repeated the calculation with the stochastic method using a second, larger value of the electromagnetic coupling $\alpha=\nicefrac{{1}}{{4\pi}}$. Using the results for the QED correction to the mass from two different values of the coupling $e$ and an ansatz $\delta m^{\textrm{stoch}}=\alpha m_{1}+\alpha^{2}m_{2}$, we can explicitly determine the $\mathcal{O}(\alpha)$- and the $\mathcal{O}(\alpha^{2})$-correction which are included in the stochastic data. More details can be found in the appendix in section C. The $\mathcal{O}(\alpha^{2})$-correction, which is included in the stochastic data, is shown by the solid green line on the right-hand side of figure 6. Thus, we find that the difference between the results from the perturbative and the stochastic method is described by the $\alpha^{2}$ contribution to the mass, which is only included in the stochastic data.

To determine the QED correction to the mass of a meson, we fit a constant to the plateau region of the QED correction to the effective mass, such as the data in figure 6. The results are given in table 1 for the stochastic and the perturbative method. We give results for charged and neutral pions as well as charged and neutral kaons. Note, that we do not include the quark-disconnected diagram for the neutral pion. One also has to keep in mind, that our calculation is not using physical quark masses, and thus the results shown here are not at the physical point. The small but significant difference in the results from the stochastic and perturbative method for charged pion and kaon is due to higher order effects in $\alpha$, which are only included in the stochastic data.

QED in a finite box is subject to substantial finite volume effects. Although in this exploratory study we do not give results at the physical point and thus, correcting for finite volume effects is not strictly necessary, we include finite volume corrections for the meson masses to illustrate that they are significant. Finite volume effects for the QED correction to the meson masses are analytically known up to $\mathcal{O}(1/L^{3})$ corrections and given by Borsanyi:2014jba

[TABLE]

with $\kappa=2.837297$. $m(L)$ and $m$ are the meson masses including QED in finite and infinite volume, respectively. In table 1 we quote results in finite volume and results $\delta m^{\textrm{inf V}}$ in infinite volume, where finite volume effects have been accounted for using equation (58).

The QED correction to the meson masses has been previously calculated in an independent calculation Blum:2010ym using a stochastic method on the same gauge ensemble, albeit different gauge configurations. A comparison of our results with the results from this independent calculation can serve as a cross check of our data. In table 2 we show results for the squared mass splitting $\Delta m^{2}=(m_{0}+\delta m)^{2}-m_{0}^{2}$ for a pion. Both light quark masses in this comparison equal the light sea-quark mass, i.e. a bare mass of $m_{l}=0.005$. We find agreement between our results and the results from this previous calculation.

4.2.2 Comparison of Statistical Errors

To compare the statistical errors on the QED correction to the effective mass between the stochastic and the perturbative data, one has to take into account, that these two datasets have not been obtained at the same numerical cost. For the stochastic data we need three inversions per quark flavour and source position ($e=0$, $e$, $-e$) to obtain the QED correction. As described in section 2.3 in our setup the calculation of the QED correction with the perturbative method requires $17$ inversions per quark flavour and source position using the single-$\mu$ insertion, if the Feynman gauge is used for the photon propagator, and $5$ inversions using the summed-$\mu$ insertion. Thus, we find a $17/3$ or $5/3$ larger numerical cost for the perturbative method to obtain the same statistics than for the stochastic data.

Figure 7 shows a comparison between the statistical errors of the stochastic and the perturbative data. The plot on the left shows the error from the perturbative data divided by the error on the stochastic data, both scaled with their respective numerical cost (i.e. the number of inversions) to have an equal cost comparison. Blue circles and purple triangles are using the single- or summed-$\mu$ insertion technique for the perturbative method, respectively. The black horizontal line shows the threshold above which the accuracy of the stochastic approach is superior to the one of the perturbative approach. We find the perturbative method to give an error which is about a factor $1.5$ to $2$ larger than the error on the stochastic method for the same costs. Comparing the two different approaches for calculating the sequential propagators for the perturbative method, we find that at the same numerical cost the statistical error is smaller when using the summed-$\mu$ insertion.

The right-hand side of figure 7 shows the ratio of the errors of perturbative and stochastic data with the same set of statistics. We find this ratio to be slightly smaller but close to one for the single-$\mu$ insertion and slightly larger then one for the summed-$\mu$ insertion, and thus finding similar statistical errors for the perturbative and stochastic data when using the same statistics for the QCD average. In summary, the ordering of statistical errors is

[TABLE]

For quenched QED, depending on whether the cost of QCD gauge configuration generation is to be included in a cost assesment, the optimal method to select will either be the most precise for the same statistics if the cost of the measurement is sub-dominant, or if a sufficient ensemble of gauge configurations already exists it makes sense to select the most precise approach for fixed measurement cost.

We note, that a cost comparison between stochastic and perturbative methods is less trivial in unquenched QED. While for the perturbative method, one needs to additionally calculate appropriate quark-disconnected diagrams, the stochastic method requires the generation of combined QCD$+$QED gauge configurations.

4.3 Strong Isospin Breaking Correction

In the following we show results for the strong isospin breaking corrections to meson masses. As discussed above, we use two different approaches to account for strong isospin breaking, one by simply using different valence up- and down-quark masses when computing valence quark propagators, and one by expanding the path integral in the quark mass difference. When comparing results from these two approaches, one has to keep in mind, that, when choosing different values for up- and down-quark masses, we fixed the up-quark mass to the isospin symmetric mass $m_{u}=\hat{m}$ and changed the down-quark mass to be $m_{d}=\hat{m}+(m_{d}-m_{u})$, where $(m_{d}-m_{u})$ approximately corresponds to the physical light quark mass difference from Fodor:2016bgu .

In the following we focus on the strong isospin correction to the masses of charged kaon $K^{+}=\overline{s}\gamma_{5}u$ and neutral kaon $K^{0}=\overline{s}\gamma_{5}d$. In this context “charged” and “neutral” refers only to the quark content, not to electromagnetic charges. In particular, we consider the strong isospin contribution to the difference $\tilde{m}_{K^{0}}-\tilde{m}_{K^{+}}$ of the masses of charged and neutral kaon. Here, we define masses denoted by $\tilde{m}$ as masses that include strong isospin corrections, but no QED effects and

[TABLE]

where $m_{m_{u}=m_{d}}$ is the isospin symmetric mass and $\delta_{\textrm{s}}m$ the strong isospin correction.

We can obtain $\tilde{m}_{K^{0}}-\tilde{m}_{K^{+}}$ by simply calculating the effective mass according to equation (49) once for a two-point function using a strange quark and a light quark with mass $m_{u}$ and once for a two-point function using a strange quark and a light quark with mass $m_{d}$ and taking their difference.

On the other hand, we obtain $\tilde{m}_{K^{0}}-\tilde{m}_{K^{+}}$ from the expansion of the path integral. According to equation (43) we find for the two-point correlation functions of charged and neutral kaon

[TABLE]

where $C_{K}^{\textrm{strongIB}}(z_{0})$ is given by (45)

[TABLE]

Since by expanding the path integral we only determine the strong isospin breaking correction which is linear in $\left(m_{d}-m_{u}\right)$, the difference $\tilde{m}_{K^{0}}-\tilde{m}_{K^{+}}$ has to be extracted using

[TABLE]

i.e. with the ratio method, as for the $\mathcal{O}(\alpha)$ corrections when using the perturbative method for QED (cf. section 4.1).

Figure 8 shows the strong isospin correction to the difference of effective masses between charged and neutral kaon. The green circles show results using different up- and down-quark masses. Here we take the difference between the cosh effective mass of a charged and a neutral kaon. The purple square points show results using the path integral expansion and equation (64).

The strong isospin breaking contribution to the difference of the masses of neutral and charged kaons is determined by fitting a constant function to the plateau region of the data shown in figure 8. From these fits we obtain

[TABLE]

for the data using different up- and down-quark masses and the path integral expansion, respectively. We find the statistical errors to be approximately the same for both methods to account for strong isospin breaking. Both methods have the same computational cost, since they require either one additional inversion with a second light quark mass or one additional inversion with a sequential insertion of the scalar current.

The authors of deDivitiis:2013xla ; deDivitiis:2011eh showed, that the strong isospin breaking correction to the mass difference between a charged and a neutral pion vanishes at $\mathcal{O}(m_{d}-m_{u})$, since the correlation functions $\tilde{C}_{\pi}^{+}=\left<\pi^{-}\pi^{+}\right>$ and $\tilde{C}_{\pi}^{0}=\left<\pi^{0}\pi^{0}\right>$ with $\pi^{+}=u\gamma_{5}\overline{d}$ and $\pi^{0}=\frac{1}{\sqrt{2}}(u\gamma_{5}\overline{u}-d\gamma_{5}\overline{d})$ receive the same leading strong isospin correction

[TABLE]

with

[TABLE]

However, the strong isospin breaking correction to the masses of neutral and charged pions differ at $\mathcal{O}\left((m_{d}-m_{u})^{2}\right)$. When calculating pion correlation functions using different input bare masses for up and down quarks, we determine strong isospin correction at all orders in $(m_{d}-m_{u})$. Using results from this approach we find a non-vanishing, albeit small mass difference

[TABLE]

5 Isospin Breaking Corrections to $a_{\mu}$

In the following, we determine the isospin breaking corrections to the anomalous magnetic moment of the muon $a_{\mu}$. In section 5.1 we give an introduction to $a_{\mu}$ and specify the setup we use to calculate the hadronic vacuum polarization. Results for the QED correction to the vector two-point function and the multiplicative renormalization $Z_{V}$ are given in sections 5.2 and 5.3, respectively. In section 5.4 we show results for the strong isospin breaking correction to $a_{\mu}$. Our results for the isospin breaking corrections to $a_{\mu}$ are summarized in section 5.5.

5.1 Introduction and Definitions

The anomalous magnetic moment of the muon $a_{\mu}$ has been experimentally measured with a precision of $\approx 0.5$ ppm at Brookhaven National Laboratory Bennett:2006fi using polarized muons in a storage ring in a magnetic field. The Standard Model estimate Olive:2016xmw of $a_{\mu}$ has been determined to the same level of accuracy. Thus, $a_{\mu}$ can serve as a high precision test of the Standard Model of particle physics. However, since many years a deviation of about $3\sigma$ persists between the experimental and theoretical estimates. This deviation might be a sign of new physics. Clearly, it is important to reduce the errors in both the experimental measurement and in the Standard Model calculation. From the experimental side, there are two upcoming experiments at Fermilab Venanzoni:2014ixa and J-PARC Otani:2015jra , both aiming to further reduce the experimental uncertainty. While the biggest contribution in the Standard Model estimate originates from the electromagnetic interaction, the largest contribution to the error comes from the strong interaction. The leading strong contribution to $a_{\mu}$ is given by the hadronic vacuum polarization (HVP).

Currently, the most precise theoretical estimate Davier:2010nc ; Hagiwara:2011af of the hadronic vacuum polarization uses the data from the cross section of $e^{+}e^{-}\rightarrow$ hadrons, and, thus, relies on experimental data. On the other hand, the HVP can be calculated from first principles using lattice QCD. The hadronic vacuum polarization $\Pi(Q^{2})$ is determined by the correlation function of two electromagnetic currents

[TABLE]

with

[TABLE]

From the HVP form factor $\Pi(Q^{2})$, the leading hadronic contribution to $a_{\mu}$ can be determined by Blum:2002ii

[TABLE]

with a kernel function $K(Q^{2})$, which is known analytically.

In recent years a lot of effort has been undertaken to determine the HVP contribution to $a_{\mu}$ using lattice calculations (see e.g Boyle:2011hu ; DellaMorte:2011aa ; Burger:2013jya ; Chakraborty:2014mwa ; Bali:2015msa ; Chakraborty:2016mwy ; Borsanyi:2016lpl ; DellaMorte:2017dyu ). However, to be competitive with the determination from $e^{+}e^{-}\rightarrow$ hadrons, a precision of $\lesssim 1\%$ is required. At this level of precision, isospin breaking corrections can no longer be neglected.

In this work we achieve the first exploratory calculation of isospin breaking corrections to the HVP, using a setup identical to the one previously described for the meson mass splittings. Note, that the QED corrections to $a_{\mu}$ are of the same order in $\alpha$ as the hadronic light-by-light scattering contribution (see Blum:2015gfa ; Blum:2016lnc ; Green:2015sra ; Asmussen:2016lse for lattice calculations of the hadronic light-by-light scattering).

For the calculation of the HVP we choose a setup with a local vector current at the source and a conserved vector current at the sink

[TABLE]

see Blum:2016xpd for further details of the framework for our calculation of the hadronic vacuum polarization. The conserved vector current $V^{c}_{\mu}$ for the Domain Wall Fermion formulation used in this work is given in equation (105). The local vector current $V^{\ell}_{\nu}$ requires a multiplicative renormalization $Z_{V}$. From the correlation function (73) we construct the HVP tensor as (see e.g. Blum:2016xpd )

[TABLE]

In (74) we have subtracted the zero-mode $\sum_{x}C_{\mu\nu}(x)$ of the vector-vector correlation function Bernecker:2011gh , which vanishes in the infinite volume limit. In Blum:2016xpd the authors showed that the zero-mode subtraction greatly reduces the statistical error on the HVP for low $Q^{2}$ when using $\mathbb{Z}_{2}$ Wall sources for the quark propagators. We also find such an improvement for the QED correction to the HVP, reducing the error for the smallest $Q^{2}$ by a factor of $\approx 4$ for the up quark and $\approx 20$ for the strange quark.

For the determination of the HVP form factor $\Pi(\hat{Q}^{2})$ we use the spatial components of (74)

[TABLE]

with vanishing spatial momentum $\vec{Q}=0$.

The QED correction to the hadronic vacuum polarization $\delta\Pi$ is determined by the QED correction to the correlation function $C_{\mu\nu}(x)$ (73). Since we use the local vector current at the source, we have to apply the appropriate multiplicative renormalization $Z_{V}$. This multiplicative renormalization $Z_{V}$ itself receives a QED correction once electromagnetism is switched on. Thus, the QED correction to the local-conserved vector two-point function $C_{\mu\nu}(x)$ is given by

[TABLE]

where $\left<V^{c}_{\mu}(x)V^{\ell}_{\nu}(0)\right>_{0}$ is the vector two-point function without QED, $Z_{V}^{0}$ the multiplicative renormalization without QED, $\delta Z_{V}$ the QED correction to $Z_{V}$ and $\delta\left<V^{c}_{\mu}(x)V^{\ell}_{\nu}(0)\right>$ the QED correction to the vector two-point function. It follows from equation (76) that the QED correction to the HVP is given by

[TABLE]

where $\delta^{Z_{V}}\Pi(\hat{Q}^{2})$ and $\delta^{V}\Pi(\hat{Q}^{2})$ are the QED corrections from the correction to $Z_{V}$ and from the correction to the vector two-point function $\left<V^{c}_{\mu}(x)V^{\ell}_{\nu}(0)\right>$, respectively. Similarly, we define the QED correction to $a_{\mu}$ as

[TABLE]

In general, $\delta a_{\mu}$ also receives a contribution from the QED correction to the lattice spacing. The lattice spacing enters in the kernel function $K(Q^{2})$ in equation (72), which depends on the muon mass. However, in this work we did not determine the lattice spacing in the presence of QED (cf. section 3).

Our results for the QED correction $\delta^{V}\!a_{\mu}$ from the vector two-point function are presented in section 5.2 and results for the QED correction $\delta^{Z_{V}}\!a_{\mu}$ from the multiplicative renormalization are given in section 5.3.

5.2 QED Correction to the Vector Two-Point Function

In this section we discuss the QED correction $\delta\left<V^{c}_{\mu}(x)V^{\ell}_{\nu}(0)\right>$ to the vector two-point function. As described above, we use a conserved vector current $V^{c}_{\mu}$ at the sink when calculating the hadronic vacuum polarization. However, the conserved current depends on the link variables $U_{\mu}(x)$ (cf. equation (105)) and thus, in the presence of QED, on the photon fields and the electromagnetic charge $e$. In the following, we refer to the conserved vector current including the $U(1)$ photon fields as $V_{\mu}^{c,e}(x)$ to indicate the dependence on $e$. Thus, for the QED correction to the HVP we have to calculate the expectation value of an operator that itself depends on the electromagnetic coupling. This has to be taken into account when expanding the path integral for the perturbative method

[TABLE]

This leads to two additional terms, that are not present in the QED correction to the meson masses. The corresponding diagrams are shown in figure 9. A detailed derivation can be found in the appendix B.2. The construction of these two terms does not require any additional inversions compared to the diagrams which we have already considered for the meson masses.

5.2.1 Results

Figures 10 and 11 show the QED correction $\delta^{V}\Pi(\hat{Q}^{2})$ to the hadronic vacuum polarization form factor for up and strange quarks, respectively. The plots on the left-hand side of both figures show results from the perturbative and the stochastic method. For the perturbative data the results shown have been calculated using the single-$\mu$ insertion, which, for the same amount of statistics, gives a smaller statistical error than the summed-$\mu$ insertion (see section 5.2.2 for a detailed comparision of statistical errors). For the multiplicative renormalization $Z_{V}^{0}$ of the local vector current we use a value determined from the ratio of the local-conserved and the local-local vector two-point functions. Further details can be found in section 5.3, where we will also determine the QED correction to $Z_{V}$.

The plots on the right-hand side of figures 10 and 11 show the correlated difference between the data from perturbative and stochastic methods. We find the data from both methods differs for a large range of $Q^{2}$ at the level of about $1-1.5\leavevmode\nobreak\ \sigma$. In order to understand this small difference we perform a computation with a second value of the electromagnetic coupling $\alpha=\nicefrac{{1}}{{4\pi}}$ for the stochastic method, so that we can distinguish between the leading and higher-order QED correction in the data from the stochastic method. We find that the deviation seen between stochastic and perturbative data for large $\hat{Q}^{2}$ to be consistent with $\mathcal{O}(\alpha^{2})$ corrections. More details on the $\mathcal{O}(\alpha^{2})$ QED corrections from the stochastic data are given in appendix C.

The calculation of $a_{\mu}$ from the HVP form factor requires the subtracted HVP $\hat{\Pi}(\hat{Q}^{2})=\Pi(\hat{Q}^{2})-\Pi(0)$. For the QED correction we therefore need to determine $\delta^{V}\Pi(0)$.

The subtracted hadronic vacuum polarization $\hat{\Pi}(\hat{Q}^{2})$ can be directly determined from the vector two-point function by Bernecker:2011gh

[TABLE]

Similarly, we calculate the QED correction to the subtracted HVP from the QED correction to the vector two-point function

[TABLE]

The results for $\delta^{V}\hat{\Pi}$ can then be used to calculate the QED correction to $a_{\mu}$ according to equation (72). We use a sine cardinal interpolation Blum:2016xpd to obtain the HVP also at non-lattice momenta

[TABLE]

where $n_{0}$ can lie anywhere in $[-T/2,T/2)$ and not only on integer values. To obtain $a_{\mu}$ we integrate using the trapezoidal rule up to momenta $\hat{Q}^{2}\approx 3$ GeV2. The integrand in the integral to obtain $a_{\mu}$ (cf. equation (72)) is peaked at small momenta around the muon mass. Contributions from momenta $>3$ GeV2 are very small, and we neglect these in this study. The results for the QED corrections to $a_{\mu}$ from the QED correction to the vector two-point function are given in table 3 alongside results without QED.

We find the QED correction $\delta^{V}\!a_{\mu}$ for the up quark to be of the order of $\lesssim 1\%$ of the value without QED. Results for the down quark can be obtained by multiplying the values for the up quark with the appropriate charge factor $1/4$ for $a_{\mu}^{0}$ and $1/16$ for $\delta^{V}a_{\mu}$. In contrast to the QED correction for the light quarks, we find the QED correction for the strange quark contribution to be negative. Although we find agreement between the HVP form factor from the perturbative and stochastic data (cf. figures 10 and 11), we find the results for the QED correction $\delta^{V}a_{\mu}$ given in table 3 to differ between the stochastic and perturbative approach by $2-3\sigma$. This is due to the $1-2\sigma$ deviation between both datasets for small $\hat{Q}^{2}$ in the HVP form factor. When calculating the subtracted HVP using equation (80) this difference gets enhanced over the whole $Q^{2}$ region for $\hat{\Pi}(\hat{Q}^{2})$.

Another method to determine $\Pi(0)$ is to fit the HVP form factor to extrapolate to $Q^{2}=0$. Suitable fit functions are given by Padé approximants Aubin:2012me

[TABLE]

In this work we use Padé $R_{11}$, which has one pole

[TABLE]

To obtain a fit function for the QED correction from (84), we allow each parameter to receive a QED correction

[TABLE]

At $\mathcal{O}(\alpha)$ we find

[TABLE]

as an ansatz for fitting the QED correction to the HVP. Since $\delta^{V}R_{11}(Q^{2})$ also depends on the parameters $a^{0}$ and $b^{0}$ from the Padé without QED, we perform a combined fit of the HVP without QED and the QED correction $\delta^{V}\Pi(\hat{Q}^{2})$. Results of these fits are shown in figure 12 for the QED correction to the HVP for the up and the strange quark. The dashed blue curve shows the fit result for the stochastic data (blue circles), the solid red curve shows the fit result for the perturbative data (red squares).

For the calculation of $a_{\mu}$ according to equation (72) we use the fit result of the Padé in the fit range, which is indicated by the range in which the Padé function is plotted in figure 12. For higher $Q^{2}$ we use the data and trapezoidal rule for the integration. As before, we integrate up to $\approx 3$ GeV2.

The results for $a_{\mu}$ without QED as well as the QED corrections from perturbative and stochastic data using the Padé $R_{11}$ are given in table 4. We find the results for the QED correction to $a_{\mu}$ for the up quark to be smaller than the values in table 3 determined using equation (80). For the strange quark we again find a negative QED correction to $a_{\mu}$, which is in agreement with the results given in table 3. We find the results from the perturbative and the stochastic data to differ by $1-2\sigma$. This difference is not as pronounced as when using equation (80) (cf. results in table 3), since the deviation between both data sets at small $\hat{Q}^{2}$ is reduced by the Padé fit as one can see on figure 12.

The difference of the results for $\delta^{V}\!a_{\mu}$ using the two different methods to determine $\hat{\Pi}(\hat{Q}^{2})$ discussed here is shown in table 5 (i.e. the difference between results from tables 3 and 4). This difference mainly arises from the large statistical errors on the QED correction for small $\hat{Q}^{2}$. A better resolution of the QED correction in this region would allow for a more reliable determination of $\Pi(0)$ and thus $\hat{\Pi}(\hat{Q}^{2})$.

Currently we do not include any quark-disconnected diagrams in the calculation of QED corrections. However, one type of quark-disconnected diagrams is also present in the electro-quenched approximation. A sketch of this diagram is shown in figure 13. Note, that in this context one is only interested in the case, where the two quark lines are additionally connected by gluons. If the two quark lines are only connected by the photon and not by gluons this diagram is conventionally counted as a higher order HVP contribution (see e.g. Jegerlehner:2009ry ), not as a QED correction to the leading order HVP. We include this diagram neither in the stochastic nor in the perturbative data.

The counterpart of this diagram without QED, i.e. without a photon coupling the two quark loops, is $SU(3)$ suppressed (see Francis:2014hoa ; Chakraborty:2015ugp ; Blum:2015you for lattice QCD calculations of the quark-disconnected contribution to the HVP and DellaMorte:2010aq ; Bijnens:2016ndo for estimates of the disconnected HVP in chiral perturbation theory). However, when including a photon to obtain the QED correction shown in figure 13 the corresponding correlation function is no longer $SU(3)$ suppressed. Thus, this quark-disconnected diagram might give a large QED correction to the HVP. We plan to include this contribution in future calculations. Note, that the diagram shown in figure 13 determines the mixing of $\rho$- and $\omega$-mesons.

5.2.2 Comparison of Statistical Errors

To compare the statistical errors between the perturbative and the stochastic data, we consider their ratio. Figure 14 shows the error on the perturbative data divided by the error on the stochastic data. For the plot on the left-hand side, the errors are scaled by the total number of inversions used in each case to obtain an equal cost comparison. Closed and open symbols denote results from the the single- or summed-$\mu$ insertion technique for the perturbative method, respectively. The horizontal black line shows “1” where both methods would give the same precision with the same numerical cost. However, we find the statistical error from the perturbative method to be larger then the error from the stochastic method. Comparing the two different approaches for calculating the sequential propagators for the perturbative method, we find that at the same numerical cost the statistical error is smaller when using the summed-$\mu$ insertion (open symbols).

The plot on the right-hand side of figure 14 shows the ratio of errors in an equal statistics comparison. We find this ratio to be smaller but close to one for the single-$\mu$ insertion and slightly larger then one for the summed-$\mu$ insertion.

We find the same ordering of statistical errors as for the QED correction to meson masses

[TABLE]

One has to keep in mind, that our study is done using unphysical quark masses and this might be a mass dependent finding. Indeed we observe a trend in an increasing ratio of errors from the perturbative over the stochastic method as the quark mass is decreased, suggesting that this ratio might even be larger for physical quark masses.

5.3 QED Correction to $Z_{V}$

The calculation of the HVP using a local current at the source (cf. equation (73)) requires the determination of the appropriate multiplicative renormalization $Z_{V}$. When including QED in the lattice calculation also $Z_{V}$ obtains an electromagnetic correction

[TABLE]

This results in a further correction $\delta^{Z_{V}}\Pi(\hat{Q}^{2})$ to the HVP at $\mathcal{O}(\alpha)$. In this work, we determine the multiplicative renormalization from local-conserved and local-local vector two-point functions. We define

[TABLE]

as the local-conserved and local-local vector two-point functions without QED and

[TABLE]

as the local-conserved and local-local vector two-point functions with QED.

The renormalization of the vector current without QED can be determined from the large time behaviour of the ratio of the local-conserved and local-local vector two-point functions

[TABLE]

This ratio is shown in figure 15 for the up and strange quark. $Z_{V}^{0}$ is determined by fitting a constant to the plateau region in the data as indicated in figure 15. The results of these fits are given in table 6.

The multiplicative renormalization of the vector current including QED is given by

[TABLE]

with the QED corrections to the local-conserved $\delta C^{lc}(t)$ and the local-local $\delta C^{ll}(t)$ vector two-point function. Equation (94) implies that the QED correction to $Z_{V}$ can be determined by

[TABLE]

The results for $\delta Z_{V}$ using equation (95) are shown in figure 16. The plot on the left shows data for the up quark, the plot on the right data for the strange quark. A constant has been fitted to the plateau region of the data to obtain $\delta Z_{V}$. The results from these fits are given in table 6 alongside the results for $Z_{V}^{0}$. We find the QED correction to $Z_{V}$ to be negative and smaller than $0.5\%$ for the up quark and even smaller for the strange quark, where the QED correction is more suppressed due to the smaller charge factor.

In table 7 we give results for the additional QED correction to $a_{\mu}$ due to the QED correction to $Z_{V}$. For the up quark we find the correction to $a_{\mu}$ from $\delta Z_{V}$ to be of the same order but with a different sign than $\delta^{V}\!a_{\mu}$, the correction from the QED correction to the vector two-point function itself (cf. results in tables 3 and 4). For the strange quark both QED corrections have the same sign, but the QED correction to $a_{\mu}$ from $\delta Z_{V}$ is about an order of magnitude bigger.

5.4 Strong Isospin Breaking Correction

To determine the strong isospin breaking corrections to $a_{\mu}$ we will in the following look at the HVP for the down quark $\Pi^{d}(\hat{Q}^{2})/q_{f}^{2}$, which is determined by the correlation function (cf. equation (73))

[TABLE]

We compute the strong isospin correction to $\Pi^{d}(\hat{Q}^{2})/q_{f}^{2}$ by either using the difference of the HVP calculated with different masses for up and down quarks, $\Pi^{m_{d}}(\hat{Q}^{2})-\Pi^{m_{u}}(\hat{Q}^{2})$, or by using the expansion of the path integral. In the latter, the strong isospin breaking correction is given by

[TABLE]

with the scalar current $\mathcal{S}$. In figure 17 the strong isospin correction to the HVP is plotted against $\hat{Q}^{2}$. Green circles show the difference of the HVP calculated using different masses for up and down quark. The purple squares show results obtained from the expansion of the path integral in the quark mass, i.e. equation (97). We find the data sets from both methods to account for strong isospin to agree with each other.

To determine the strong isospin breaking correction $\delta_{\textrm{s}}a_{\mu}$ to the anomalous magnetic moment of the muon from the data shown in figure 17, we use either equation (80) to determine $\hat{\Pi}(\hat{Q}^{2})$ or a Padé fit. The Padé fit can be done in a similar way as for the QED correction, i.e. assigning a strong isospin breaking correction to each of the parameters in the Padé function. Thus, for $R_{11}$ we obtain

[TABLE]

as an ansatz to fit the data for the strong isospin correction to the HVP. The results of these fits are shown in figure 17 by the solid green line for the data using different masses for up and down quark and the dashed purple line for the data using the path integral expansion.

The results for the strong isospin breaking correction to $a_{\mu}$ are given in table 8. Using the Padé fits we are able to resolve the strong isospin breaking correction to $a_{\mu}$. When comparing the results in table 8 with the value in the isospin symmetric limit (cf. table 4) $a^{u}_{\mu}/q_{u}^{2}=(716\pm 25)\times 10^{-10}$, we find that the strong isospin correction is $\delta_{\textrm{s}}a_{\mu}/a^{u}_{\mu}\approx 0.9\%$.

In addition, we have to determine the strong isospin breaking correction to the multiplicative renormalization $Z_{V}$, which can be obtained by comparing results for $Z_{V}$ using either $m_{u}$ or $m_{d}$ as the valance quark mass. We find this correction to be very small $\delta_{\textrm{s}}Z_{V}/Z_{V}=(Z^{m_{u}}_{V}-Z^{m_{d}}_{V})/Z_{V}\approx 0.02\%$.

5.5 Summary IB Corrections to $a_{\mu}$

Our results for the QED corrections to the anomalous magnetic moment of the muon are summarized in table 9. The first column $a_{\mu}^{0}$ shows the result in the isospin symmetric limit for up and strange quarks (Note, that in the isospin symmetric limit, the contribution from the down quark is simply $1/4$ of the up quark). Results for the QED correction $\delta^{V}a_{\mu}$ from the vector two-point function are given in columns two to five and the QED correction $\delta^{Z_{V}}a_{\mu}$ from the multiplicative renormalization in columns six and seven. For the QED correction $\delta^{V}a_{\mu}$ we have determined the subtracted vacuum polarization $\hat{\Pi}(\hat{Q}^{2})$ with two different methods, either using equation (80) or using Padé $R_{11}$ to obtain $\Pi(0)$ and we quote both results separately. We find results from both techniques for determining $\hat{\Pi}(\hat{Q}^{2})$ to differ especially for the up quark. This difference mainly arises from the large statistical errors on the QED correction for small $\hat{Q}^{2}$. A reduction of the statistical error in the low $\hat{Q}^{2}$ region is required to achieve a more reliable determination of $\Pi(0)$.

The total QED correction to $a_{\mu}$ is given by the sum of the two contributions $\delta^{V}a_{\mu}$ and $\delta^{Z_{V}}a_{\mu}$. We have not added these contributions in order to illustrate, that the statistical error is dominated by the QED correction $\delta^{V}a_{\mu}$ that originates from the QED correction to the vector two-point function, while the QED correction to the multiplicative renormalization $Z_{V}$ is determined very precisely.

We find the overall QED correction to $a_{\mu}$ to be smaller than $1\%$ for the up quark, where the QED contribution is enhanced by the charge factor compared to the down and strange quark. For the strange quark we find the QED correction to be about $0.1\%$ of the isospin symmetric result.

The above findings are for unphysical sea and valence light quark masses and QED corrections to $a_{\mu}$ might be larger at the physical point.

Our results for the strong isospin breaking correction are summarized in table 10. We have accounted for strong isospin breaking by either using different masses for the valence up and down quark or by using a path integral expansion in $(m_{u}-m_{d})$. For both datasets we have determined the subtracted vacuum polarization $\hat{\Pi}(\hat{Q}^{2})$ with two different methods, either using equation (80) or using Padé $R_{11}$ to obtain $\Pi(0)$ and both results are quoted separately in table 10. We find the strong isospin correction to be $0.9\%$ of the isospin symmetric result.

6 Conclusions and Outlook

In this work we have calculated the isospin breaking corrections to meson masses and the hadronic vacuum polarization in an exploratory study on a $64\times 24^{3}$ lattice with an inverse lattice spacing of $a^{-1}=1.78$ GeV using unphysical quark masses.

We have included electromagnetic effects in the lattice calculation with two different approaches, a stochastic and a perturbative approach. To our knowledge, this work is the first direct comparison of results obtained from these two methods. In both methods, we have treated QED in an electro-quenched setup, i.e. we have considered the sea quarks as electrically neutral.

As a starting point for comparing the stochastic and perturbative methods we have calculated the QED correction to meson masses. We have shown, that these QED corrections have to be extracted differently from stochastic or perturbative data, taking into account, that the stochastic data contains QED corrections to the correlation functions from all orders in $\alpha$, while the perturbative data only include $\mathcal{O}(\alpha)$ corrections. We find the results from the perturbative and the stochastic method to be consistent with each other up to small deviations, which are of $\mathcal{O}(\alpha^{2})$. Albeit the $\mathcal{O}(\alpha^{2})$ corrections are small, we are able to resolve these with the statistics used in this study.

In this work we have determined for the first time the QED corrections to the hadronic vacuum polarization and its contribution to the anomalous magnetic moment of the muon $a_{\mu}$. We have calculated the QED correction to the local-conserved vector two-point function and to the multiplicative renormalization $Z_{V}$ for the local vector current used in our setup to calculate the HVP. An overview over our results for the QED correction to the HVP is presented in table 9. In total, we find the QED correction to $a_{\mu}$ to be $<1\%$ for the up quark and $0.1\%$ for the strange quark. However, one has to keep in mind that this calculation has not been done using physical quark masses. In addition, we have determined the strong isospin correction to the HVP, which we find to be $\approx 0.9\%$. An important conclusion from this is, that when aiming at a calculation of $a_{\mu}$ with a precision of $1\%$, QED and strong isospin breaking corrections would need to be included.

Our data allows us to directly compare the statistical precision obtained from the stochastic and the perturbative method. We find that for the QED correction to the meson masses as well as for the HVP the stochastic method results in a statistical error which is about a factor of $1.5-2$ smaller than the statistical error from the perturbative method for the same numerical cost. Thus, the stochastic method is favourable for the particular choice of simulation parameters and quantities considered in this work. However, this might differ for a study using unquenched QED, where a cost comparison between both methods is less trivial.

In this work we have not made any attempt to include finite volume corrections for the QED correction to the HVP. Finite volume corrections with photons in a finite box can be substantial (cf. e.g. the results in Borsanyi:2014jba and in table 1) and thus need to be taken into account. We are currently investigating the finite volume corrections to the QED correction to $a_{\mu}$ to include those in our calculations.

Having successfully completed this exploratory study with unphysical quark masses, a calculation of the QED corrections to the HVP at physical quark masses is under active investigation. This present work has demonstrated the methods and feasibility to enable the future work.

Acknowledgements

The authors warmly thank F. Sanfilippo for usefull discussions. This work has received funding from the STFC Grant ST/L000296/1, the EPSRC Centre for Doctoral Training in Next Generation Computational Modelling grant EP/L015382/1 and from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement 279757. V.G. acknowledges partial support from the Horizon 2020 INVISIBLESPlus (H2020-MSCA-RISE-2015 -690575). P.A.B. and A.P. are supported in part by UK STFC grant ST/L000458/1. C.L. is supported in part by US DOE Contract $\#$AC-02-98CH10886(BNL) and in part through a DOE Office of Science Early Career Award. This work used the DiRAC Blue Gene Q Shared Petaflop system at the University of Edinburgh, operated by the Edinburgh Parallel Computing Centre on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K000411/1, STFC capital grant ST/H008845/1, and STFC DiRAC Operations grants ST/K005804/1 and ST/K005790/1. DiRAC is part of the National E-Infrastructure. Computing support for this work came partially from the Lawrence Livermore National Laboratory (LLNL) Institutional Computing Grand Challenge program.

Appendix A Domain Wall Action and Currents

The Domain Wall fermion action used in this work is given by Furman:1994ky ; Blum:2000kn

[TABLE]

with

[TABLE]

and

[TABLE]

To include couplings of photon fields to the quarks in the fermionic action (99) one has to replace the gauge links in (102) and (103) as follows

[TABLE]

The conserved vector current $V_{\mu}^{c}(x)$ and the tadpole operator $T_{\mu}(x)$ are given by

[TABLE]

[TABLE]

Appendix B Expansion of the Path Integral

The expansion of the expectation value of an observable $O$ in the electromagnetic coupling $e^{2}$ is given as

[TABLE]

The leading order electromagnetic correction is thus determined by

[TABLE]

In the electro-quenched approximation we do not include QED in the fermion determinant $\det D[U,A]\equiv\det D_{0}[U]$, and consequently the partition function $Z$ does not depend on the electromagnetic coupling $e$.

B.1 Meson Two-Point Functions

For a meson two-point function the observable $O$ is of the form $O=(q_{f}\Gamma q_{f}^{\prime})(q_{f}\Gamma^{\prime}q_{f}^{\prime})$ for two quark flavours $f$ and $f^{\prime}$, and does not depend on the elementary charge $e$. Thus, equation (108) can be written as

[TABLE]

In the following we drop the dependence of $S_{F}$ on the fields $\Psi,\overline{\Psi},A,U$ for simplicity, i.e. $S_{F}\equiv S_{F}[\Psi,\overline{\Psi},A,U]$. The derivative in (109) can be written as

[TABLE]

For the Domain Wall action (99) including QED (cf. (104)) one finds

[TABLE]

and

[TABLE]

Inserting this into (109) yields

[TABLE]

with the conserved vector current (105) and the tadpole operator (106). Using

[TABLE]

one finds for the expansion (107) of the path integral at $\mathcal{O}(\alpha)$

[TABLE]

where $\left<\cdot\right>_{0}$ is the expectation value over fermionic and gluonic fields.

B.2 HVP

For the QED correction to the HVP in the perturbative method we have to expand the path integral for an operator of the form

[TABLE]

with a local vector current $V_{\nu}^{l}$ and a conserved vector current $V_{\mu}^{c,e}$, which, including QED, is given by

[TABLE]

Taking into account the explicit dependence of the operator on the electromagnetic coupling, one has to calculate

[TABLE]

The derivatives of the conserved vector current $V_{\mu}^{c,e}$ with respect to $e$ are given by

[TABLE]

and

[TABLE]

Thus, in total, one finds for the expansion of the path integral for the operator (116)

[TABLE]

The second term on the right-hand side of the first line of (121) is the tadpole diagram, the term on the second line gives rise to the quark self-energy and the photon exchange diagram. The terms in the third and fourth line are the two terms shown in figure 9 and originate from the expansion of the operator.

B.3 Strong Isospin Breaking

For determining the strong isospin correction using the path integral expansion, we have to calculate

[TABLE]

The derivative of the expectation value with respect to $e$ is given by

[TABLE]

For the Domain Wall Fermions used in this work, we find

[TABLE]

with four dimensional fields $\overline{\psi}$, $\psi$. Thus, we find for the path integral expansion (122)

[TABLE]

with the scalar current $\mathcal{S}=\sum\limits_{x}\overline{\psi}(x)\psi(x)$.

Appendix C $\mathcal{O}(\alpha^{2})$ Effects

C.1 Meson Masses

The results from the stochastic and the perturbative method are expected to differ by effects which are of $\mathcal{O}(\alpha^{2})$. Indeed, when comparing the QED correction to the effective mass (cf. figure 6), we find a deviation between both datasets of about $1\%$ of the QED correction itself. This deviation is consistent with $\mathcal{O}(\alpha^{2})$ corrections. To explicitly check this, we have repeated the calculation with the stochastic method for a second larger value of the electromagnetic coupling $\alpha=1/4\pi$.

Figure 18 shows the QED correction to the mass of the charged kaon for the stochastic data calculated with both values of $\alpha$ plotted against $e^{2}=4\pi\alpha$. The solid green curve is a quadratic function of the form $a\cdot e^{4}+b\cdot e^{2}$, that was matched to the two data points. The dotted grey line is the linear term $b\cdot e^{2}$ of the quadratic curve, i.e the leading order QED contribution. The plot on the right-hand side of figure 18 shows a zoom around the physical value of the coupling. One can clearly see the difference of the leading order contribution $b\cdot e^{2}$ and the full results from the stochastic method, which we find at physical $e^{2}$ to be at the same order as the statistical error. Thus, with the statistics used in this work, we are able to resolve also the $\alpha^{2}$ effects.

The $\mathcal{O}(\alpha^{2})$ effect that we obtain from the stochastic data is shown by the solid green line on the right-hand side of in figure 6. We find the $\mathcal{O}(\alpha^{2})$ contribution to be consistent with the difference that we observe between stochastic and perturbative data.

In Matzelle:2017qsw the authors discuss finite volume effects for the next-to-leading order QED corrections to meson masses. However, for our calculation with the physical value of $\alpha$ we find the $\mathcal{O}(\alpha^{2})$ effects to be about $1\%$ of the leading QED correction, and thus, considering different finite volume effects for the $\mathcal{O}(\alpha^{2})$ contributions included in the stochastic data is not relevant at the current level of precision.

C.2 HVP

When comparing the results for the HVP from perturbative and stochastic data (cf. figures 10 and 11) we found a deviation between both datasets at the level of $1-1.5\sigma$. To check if this deviation originates from $\mathcal{O}(\alpha^{2})$ effects, which are only included in the stochastic data, we use results from a computation with a larger value of the coupling $\alpha=1/4\pi$. Using the data from two different values of $\alpha$ for the stochastic method, we are able to extract the $\mathcal{O}(\alpha)$ contribution at the physical value of the coupling. For every value of the four-momentum transfer $Q^{2}$ we match a quadratic curve $a\cdot e^{4}+b\cdot e^{2}$ through the two data points, similarly as described for the QED correction to the meson masses above. The leading order QED correction is determined by the term linear in $e^{2}$. The results for the leading QED correction from the stochastic data can then be compared with the results from the perturbative method. Figure 19 shows the correlated difference between stochastic and perturbative data for up quarks (left plot) and strange quarks (right plot). Purple circles show the difference using the results from the stochastic data which still include effects to all orders in $\alpha$ (i.e. the same points that were already shown in figures 10 and 11). The light blue triangles show the difference using the $\mathcal{O}(\alpha)$ contribution from the stochastic data. We find the $\mathcal{O}(\alpha)$ data from the stochastic method to be in agreement with the results from the perturbative method over a wide range of $Q^{2}$. Thus, the difference between stochastic and perturbative data found in section 5.2.1 is consistent with effects which are of higher order in $\alpha$.

When calculating the QED correction to $a_{\mu}$ using only the $\mathcal{O}(\alpha)$ contribution from the stochastic data, we find the change in $\delta^{V}\!a_{\mu}$ to be much smaller than the statistical errors itself. Thus, we find $\mathcal{O}(\alpha^{2})$ effects to be not relevant for $\delta a_{\mu}$ at this level of precision.

Appendix D Comparison Coulomb and Feynman Gauge

The quantities calculated in this work (QED corrections to meson masses and HVP) are expected to be gauge invariant. Gauge invariance is not broken by in the QEDL prescription, i.e. subtracting the spatial zero modes, since gauge invariance is realized separately on every mode in momentum space. To numerically check for gauge invariance, we compare results using the Feynman gauge with results using the Coulomb gauge. Coulomb gauge photon fields can be obtained from Feynman gauge photon fields using an appropriate projector Borsanyi:2014jba

[TABLE]

The photon propagator in Coulomb gauge is given by

[TABLE]

with $\tilde{k}_{\mu}\equiv(0,\vec{k})_{\mu}$. The phase factor $\exp(ik\cdot(\hat{\mu}-\hat{\nu})/2)$ in equation (127) originates from the Fourier transformation with photon fields defined on the mid-links of the lattice. Note that this phase factor cancels for diagonal contributions $\mu=\nu$.

For the stochastic method we have calculated the QED contributions using the same set of statistics with the Feynman and the Coulomb gauge. For the perturbative method the calculation using the Coulomb gauge is more expensive in our setup, since also contributions from $\mu\neq\nu$ have to be determined. Thus, we restrict the calculation using the Coulomb gauge with the perturbative method to only one source position.

In figure 20 the QED correction to the effective mass of a charged kaon is shown using the Coulomb and the Feynman gauge for the stochastic method (left) and the perturbative method (right). Purple triangles show results in the Feynman gauge, orange circles results in the Coulomb gauge.

We find agreement between the QED correction to the effective mass in the Feynman and the Coulomb gauge for large $t$, where the QED correction to the meson mass is determined. For small $t$ where the data contains contributions from excited states, we find deviations between Feynman and Coulomb gauge. The data in this region also depends on the creation amplitudes of the states, which are not necessarily gauge independent quantities.