Lattice QCD Calculation of the Pion Mass Splitting

We use the infinite volume reconstruction method to calculate the charged/neutral pion mass difference. The hadronic tensor is calculated on lattice QCD and then combined with an analytic photon propagator, and the mass splitting is calculated with exponentially-suppressed finite volume errors. The calculation is performed using six gauge ensembles generated with $2+1$-flavor domain wall fermions, and five ensembles are at the physical pion mass. Both Feynman and Coulomb gauge are adopted in the calculation and result in a good agreement when the lattice spacing approaches zero. After performing the continuum extrapolation and examining the residual finite-volume effects, we obtain the pion mass splitting $\Delta m_\pi = 4.534(42)(43)~\mathrm{MeV}$, which agrees well with experimental measurements.


INTRODUCTION
One of the central goals in high energy physics is understanding the nature of the matter that we observe in the universe. As one of the four known fundamental interactions, the strong interaction binds together quarks and gluons into hadrons and most of the hadron mass arises in turn from the binding energy; the individual quarks provide only a very small portion of the mass. Precision calculation based on the lattice formulation of quantum chromodynamics (QCD), the theory of the strong interaction, can manifest its success by computing the various hadron spectra, which agree well with experimental measurements [1]. When the precision reaches percent or subpercent level, another fundamental force, electromagnetic interaction, is urged to be considered in theoretical calculations, although its effects are suppressed by a factor of the fine-structure constant α EM ≈ 1 137. Inclusion of the electromagnetic corrections does not only provide precise hadron spectroscopy [2][3][4][5][6][7][8][9][10][11][12], but also plays an important role in studies of leptonic and semileptonic decays of hadrons [13][14][15][16][17][18][19], which dramatically expands the horizon of lattice QCD studies.
Among various hadrons, pions play a unique role in the development of theoretical particle and nuclear physics. Pions were first proposed by H. Yukawa in 1935 as the carrier particles of the strong nuclear force [20]. As Nambu-Goldstone bosons [21,22], pions result from the spontaneous break down of chiral symmetries of QCD effected by quark condensation and serve as the active degrees of freedom sensitive to chiral dynamics [23][24][25]. Besides, the anomalous decay rate of the neutral pion led to the discovery of Adler-Bell-Jackiw anomaly of quantum electrodynamics (QED) [26,27], which revealed for the first time the violation of classical symmetry by quantum corrections. It can be concluded that the thorough study of the nature of the pions is a key to our better un-derstanding of QCD and the strong interaction. In this work, we focus on the study of the mass splitting between the charged and neutral pions, which represents the interplay between two fundamental interactions, the strong and the electromagnetic. For two reasons, the pion mass splitting is ideal for a lattice QCD calculation and for an exploratory study of new methodology. First, pions are the lightest hadrons and their correlation functions have very good statistical signals. Second, the isospin breaking effects of the up and down quark mass difference are suppressed by a factor of (m u − m d ) 2 Λ 2 QCD ∼ 10 −4 with m u d the up/down quark mass and Λ QCD ≈ 350 MeV the nonperturbative QCD scale, leaving the electromagnetic effects the leading contribution to the pion mass splitting. Thus the ambuguity in separating the isospin breaking effects from the electromagnetic interaction and the quark mass difference becomes irrelevant in this study. One can simply perform the lattice QCD calculation in the isospin symmetric limit and compute the four-point correlation functions for the QED self-energy diagrams.
In practice we adopt the infinite-volume reconstruction (IVR) method proposed in Ref. [28], which allows us to calculate electromagnetic corrections to stable hadron masses with only exponentially suppressed finite-volume effects. Upon its development, this method has been successfully applied and extended to various electroweak processes involving photon or massless leptonic propagators [18,[29][30][31][32]. This is the first time that we apply this methodology to the calculation of the pion mass splitting. The calculation includes the complete diagrams with both connected and disconnected quark-field contractions. We employ two gauge fixings for the photon propagator and confirm that the lattice results are consistent in the continuum limit. We utilize the random field sparsening technique proposed in Ref. [33], which allows us to improve the precision of the correlation functions significantly with only a modest cost of computational Reference m π ± − m π 0 (MeV) RM123 2013 [5] 5.33(48)stat(59)sys a R. Horsley et al. 2016 [7] 4.60(20)stat RM123 2017 [9] 4.21 (23)  resources. By using five gauge ensembles generated with N f = 2 + 1 domain wall fermions at physical pion mass and one additional ensemble at m π ≈ 340 MeV, we obtain the pion mass splitting with a percent-level uncertainty, which is about 5-10 times smaller than previous lattice QCD calculations [5,7,9,34,35]. (See Table I. Refs. [34,35] present the pioneering quenched calculations and thus the results are not included in Table I.) For the first time in the literature, we have clearly resolved and included the contribution from the quark disconnected diagram to the pion mass splitting (see Eq. (17) and the diagram below this equation). This diagram is related to the π 0 − η − η ′ mixing and has also been calculated in Ref. [36].

INFINITE VOLUME RECONSTRUCTION METHOD
The hadron mass extraction relies on the calculation of the hadron QED self-energy for a stable hadronic state N via the following infinite volume Euclidean space-time integral where the hadronic part H µ,ν (x) = H µ,ν (t, ⃗ x) is given by where N (⃗ p)⟩ indicates a hadronic state N with mass M and spatial momentum ⃗ p, J µ = 2eūγ µ u 3 − edγ µ d 3 − esγ µ s 3 is the electromagnetic current, and S γ µ,ν is the photon propagator whose form is analytically known. In Ref. [28], we introduced the IVR method to relate the infinite volume integration of infinite volume hadronic matrix elements H µ,ν to a finite lattice volume integration of finite volume matrix elements H L µ,ν with only exponentially suppressed finite volume errors. This is accomplished with the following three steps: 1. We pick t s to separate the infinite volume integral into two parts: where I (s) and I (l) are the "short distance" ( x t < t s ) and "long distance" ( x t ≥ t s ) contributions respectively.
2. For sufficiently large t s , I (l) is dominated by the lightest single particle intermediate states, and can be calculated using the hadronic matrix elements at fixed time separation t s . The excited-state effects ignored in this step are exponentially suppressed by large t s .
3. The next step is to approximate H µ,ν using H L µ,ν and restrict the integration region to the finite lattice volume. This step only introduces exponentially suppressed finite volume errors since we only need to use H L µ,ν (x) for x t ≤ t s ≲ L. The final formula obtained in Ref. [28] is expressed in terms of lattice-calculable quantaties: where L µν (t s , ⃗ x) is an infinite volume QED weighting function, defined as: In this paper, we use t s = L 2 for our final result to ensure both the excited states and finite volume errors introduced in steps 2 and 3 described above are exponentially suppressed as we increase the lattice size L.

Gauge-Specific Expressions
We also present the relevant expressions for the photon propagator S γ µ,ν (x) and QED weighting function L µ,ν (t s , ⃗ x) in Feynman and Coulomb Gauges: At the origin (x = 0), the continuum photon propagator is divergent. In our lattice calculation, we regularize this divergence with the following choice:

RELEVANT DIAGRAMS
Contributions to the charged/neutral QED pion mass splitting are derived from the following hadronic matrix element: where π( ⃗ 0) represents a pion with zero momentum. Only two contractions contribute to this matrix element [5]. One yields (where S represents the light quark propagator) C 1 µ,ν (x − y) = ⟨Tr S(x; t src )γ 5 S(t src ; x)γ µ (17) ×Tr S(y; t snk )γ 5 S(t snk ; y)γ ν ⟩ QCD which is related to the following quark disconnected diagram.

γ π π
The other possible contraction yields (up to the insertion of a photon propagator) which is represented by the following quark connected diagram. γ π π To ensure projection onto the pion state, we fix the time separation between the pion interpolating operators and the closest electromagnetic current operators to be t sep = t snk − x t = y t − t src for both diagrams (assuming x t ≥ y t ). The values of t sep for each ensemble are listed in Table II. Coulomb gauge fixed wall source is used for the pion interpolating operators. Combining these diagrams yields the hadronic contribution to the mass shift, which is represented by C π (t 2 − t 1 ) = ⟨Tr S(t 2 ; t 1 ))γ 5 S(t 1 ; t 2 )γ 5 ⟩ QCD .
Local lattice vector current operators are used in the contraction, and the corresponding renormalization factor Z V is included in the above formula. We have taken the around the world effects into account when calculating the pion correlation function in Eq. (20). In principle, both the isospin breaking due to the up / down quark mass difference and the QED effects described above contribute to the pion mass difference. For many hadronic observables, the corrections due to the up / down quark mass difference are at the same order as the QED corrections. Therefore, we treat O((m d − m u ) Λ QCD ) corrections to be of the same perturbative order as corrections of order O(α QED ). For pion mass splitting, the O((m d −m u ) Λ QCD ) effect is zero. [5] Therefore, to order O(α QED , (m d − m u ) Λ QCD ), the two diagrams discussed above represent the only contributions to the pion mass splitting.

NUMERICAL RESULTS
Calculation of the hadronic function H L µν (x) was performed on six ensembles generated by the RBC and UKQCD collaborations. The names and attributes of the ensembles are shown in Table II. The properties of these ensembles are studied in detail in Ref. [37]. In particular, the lattice spacing and other basic parameters of the ensemble are determined by matching the lattice calculation of the masses of pion, kaon, and the Ω baryon to their experimental values. We use All-modesaveraging (AMA) [38,39], zMöbius [40], and compressed eigenvector deflation [41] methods to accelerate the calculation of the propagators. Figure 1 shows the mass shift ∆m π ≡ m π + − m π 0 as a function of t s in the Feynman and Coulomb gauges for the 24D and 32D ensembles. It can be seen from the plots that the finite volume effects are very small as the differences between 24D and 32D are barely visible. Also, for t s ≳ 1.5 fm, the final results have a very mild dependence on t s , suggesting the excited states contribution beyond t s , which is exponentially suppressed and ignored in the IVR method, is indeed quite small. In the following analysis, we stick to: With this choice, the finite volume effects at fixed t s and the excited intermediate states' contribution beyond t s are both exponentially suppressed by the spatial lattice size L, and we will refer to the sum of these two effects as the finite volume effects in the following analysis. In Feynman gauge, the difference between 32D and 24D is −0.035 (16) MeV. This is consistent with a scalar QED calculation, which yields −0.022 MeV [28]. In Coulomb gauge, the difference between 32D and 24D is 0.002 (17) MeV.
The results for each ensemble are presented in Table III. In the table, we also show the Coulomb potential contribution to the pion mass difference. This contribution comes from the time component of the Coulomb gauge photon propagator, S γ,Coul t,t in Eq. (9).

EXTRAPOLATION TO THE PHYSICAL POINT
We perform the continuum extrapolation for the Iwasaki ensembles (48I and 64I) and the I-DSDR ensembles 24D, 32Dfine, and 24DH separately using the following formula (m π,phys = 135 MeV). m π ∆m π (a 2 , m π ) = m π,phys ∆m π (0, m π,phys ) (22) When using Domain wall fermions in lattice calculations, lattice artifacts, which scale as O(a) or O(a 2n+1 ), are absent [37,42]. We choose to include an O(a 2 ) term to account for the lattice artifact. We also include a naive pion mass dependence term in this extrapolation formula to accommodate the small pion mass difference in the I-DSDR ensembles. This pion mass dependence term does not effect our main result, which is taken from the continuum extrapolation of the Iwasaki ensembles exactly at the physical pion mass. The extrapolation results of the Iwasaki 48I and 64I ensembles are shown in Figure 2, where separate fits were conducted for Feynman and Coulomb gauges. To estimate the discretization systematic error for the fitted value of ∆m π (0, m π,phys ), we perform two slightly different fits by: 1. Using the same fitting formula Eq. (22), but excluding the contribution from x = 0 for all ensembles.
2. Using a different fitting formula which includes an additional O(a 4 ) term: m π ∆m π (a 2 , m π ) = m π,phys ∆m π (0, m π,phys ) (23) where the magnitude of the O(a 4 ) term is assumed to equal to the square of the O(a 2 ) term as an estimate of the remaining systematic effects.
We obtain the differences between the results of the above fits and the original fits. The maximum of the two differences is used as the estimation of the remaining discretization systematic error.
After continuum extrapolation, we use the differences between 32D and 24D to correct the finite volume effects. The absolute size of the correction is used as the estimation of the remaining finite volume systematic error. The continuum extrapolated, finite volume corrected results are shown in Table IV and Table V. The discretization and finite volume systematic errors are combined in quadrature. As expected, the continuum extrapolations from the I-DSDR ensembles (Table V) have larger discretization systematic errors due to larger lattice spacings. Therefore, we use the continuum extrapolation from the Iwasaki ensembles (Table IV) as our main results.
As a byproduct of the calculation, we plot the Coulomb potential contribution as a function of spatial separation  at t = 0 in Figure 3. This plot provides some indication of the size and the shape of the pion. CONCLUSION We have used the IVR method to calculate the pion mass splitting to order O(α QED , (m u − m d ) Λ QCD ) on the lattice. The calculation is directly performed at the physical pion mass (m π = 135 MeV), with different lattice spacings and lattice sizes. Both the connected and  The Coulomb potential contribution to the pion mass difference. The curve is the partial sum, I C (x) = 1 2 ∫ dt ∫ ⃗ y ≤x d 3 ⃗ y H L t,t (t, ⃗ y)S γ t,t (t, ⃗ y). Error bars are statistical only. We use the results from the 48I and 64I ensembles to obtain the continuum limit, and include the finite volume corrections from the 32D and 24D difference. the disconnected diagrams are included. We have obtained the continuum and infinite volume limit results with Feynman gauge: ∆m π = 4.534(42)(43) MeV, where the quantities in the first and second sets of parenthesis correspond respectively to the statistical and systematic uncertainties. We have also performed the same calculation using Coulomb gauge. Results in both gauges are consistent with each other and with the experimental value 4.5936(5) MeV [43]. Such high precision agreement demonstrates the success of both lattice QCD and the IVR methods. We plan to use the IVR method in future lattice QCD + QED spectroscopy calculations. DOE grant DE-SC0010339. X.F. Working with the second line of Eq. 17 and taking the leftmost partial derivative first gives: Taking the second partial derivative shows that: Putting it all together and restoring the (µ, ν) = (t, t) contribution we find