# Strong-isospin-breaking correction to the muon anomalous magnetic moment from lattice QCD at the physical point

###### Abstract

All lattice-QCD calculations of the hadronic-vacuum-polarization contribution to the muon’s anomalous magnetic moment to date have been performed with degenerate up- and down-quark masses. Here we calculate directly the strong-isospin-breaking correction to for the first time with physical values of and and dynamical , , , and quarks, thereby removing this important source of systematic uncertainty. We obtain a relative shift to be applied to lattice-QCD results obtained with degenerate light-quark masses of = +1.5(7)%, in agreement with estimates from phenomenology.

###### pacs:

^{†}

^{†}preprint: FERMILAB-PUB-17-486-T

Fermilab Lattice, HPQCD, and MILC Collaborations

## I Introduction

The anomalous magnetic moment of the muon, , provides a stringent test of the Standard Model and a sensitive probe of new particles and forces beyond. It has been measured by BNL Experiment E821 to a precision of 0.54 parts-per-million Bennett et al. (2006), and the experimental result disagrees with Standard-Model theory expectations by more than three standard deviations Blum et al. (2013). To investigate this discrepancy, the Muon Experiment at Fermilab aims to reduce the experimental error by a factor of four, with a first result competitive with E821 expected within a year Grange et al. (2015). To identify definitively whether any deviation observed between theory and experiment is due to new particles or forces, the theory error must be reduced to a commensurate level.

The largest source of theory uncertainty is from the O() hadronic vacuum-polarization contribution Blum et al. (2013), which is shown in Fig. 1 and is denoted by throughout this work. A well-established method for determining this contribution employs dispersion relations combined with experimentally measured electron-positron inclusive scattering cross-section data. Recent determinations from this approach quote errors of 0.4–0.6% Jegerlehner (2017a); Keshavarzi et al. (2017); Davier et al. (2017). Numerical lattice quantum chromodynamics (QCD) provides a complementary, systematic approach for calculating directly from first-principles QCD. Several independent lattice-QCD efforts to obtain are ongoing Aubin and Blum (2007); Chakraborty et al. (2017); Della Morte et al. (2017); Lellouch (2017); Lehner (2017), with errors on recent determinations ranging from about 2–6% Burger et al. (2014); Chakraborty et al. (2017); Della Morte et al. (2017); Lellouch (2017). The theoretical precision on needed to match the target experimental uncertainty is about 0.2%.

In this Letter, we remove one of the largest systematic errors common to all current lattice-QCD calculations of , namely that arising from the use of degenerate up- and down-quark masses. We do so by calculating directly the strong-isospin-breaking contribution to the light-quark-connected contribution to . Phenomenological estimates suggest that the effect of strong-isospin breaking on is about 1% Wolfe and Maltman (2011); Jegerlehner and Szafron (2011); Jegerlehner (2017b). Electromagnetic effects are also not yet included in lattice-QCD calculations of and lead to a similar-sized uncertainty Hagiwara et al. (2004); Davier et al. (2010). In order to disentangle quark-mass from electromagnetic effects, we define the strong-isospin-breaking correction using up- and down-quark masses tuned from experimental hadron masses with QED effects removed Basak et al. (2016).

The effect of strong-isospin breaking on the light- and strange-quark connected contributions to has been calculated in an exploratory study by the RBC/UKQCD Collaborations Boyle et al. (2017) in three-flavor lattice-QCD, with a heavy pion mass of about 340 MeV, and isospin-symmetric sea quarks. Preliminary four-flavor results for the strong-isospin-breaking contribution to have also been presented by the ETM Collaboration Giusti et al. (2017) for several pion masses down to MeV, but with low statistics. In this Letter, we analyze two QCD gauge-field ensembles recently generated by the MILC Collaboration with four flavors of highly improved staggered (HISQ) sea quarks and very high statistics; see Ref. Bazavov et al. (2013) for methodology. One of the ensembles has fully nondegenerate quark masses with the and quarks all fixed to their physical values. Our calculation is the first determination of at the physical pion mass and with sea isospin breaking.

This Letter is organized as follows. We first present our numerical lattice-QCD calculation in Sec. II. Next, in Sec. III, we calculate the strong isospin-breaking correction to and discuss the contributions to the systematic error. We present our final result and compare it with phenomenological estimates and previous lattice-QCD calculations in Sec. IV.

## Ii Lattice calculation

We calculate the strong-isospin-breaking correction to on two new QCD gauge-field ensembles generated by the MILC Collaboration with four flavors of HISQ quarks Follana et al. (2007); Bazavov et al. (2013). Table 1 summarizes key parameters of the configurations. The two ensembles have the same lattice spacing, which is approximately 0.15 fm, and the same strange- and charm-quark masses, which are both fixed close to their physical values. With staggered quarks, the pions possess an additional “taste” quantum number. Discretization errors from the HISQ action generate corrections to the squared sea-pion masses of different tastes. On both ensembles, the mass of the taste-Goldstone pion is fixed close to Nature’s value of MeV, which is the mass that the charged pion would have in the absence of electromagnetism. The root-mean-squared pion mass (averaging over tastes) is about MeV.

(GeV) | |||
---|---|---|---|

1.13215(35) | 0.2426/0.2426/6.73/84.47 | 1902 | 0.1347(7) |

1.13259(38) | 0.1524/0.3328/6.73/84.47 | 4963 | 0.1346(7) |

The two ensembles differ in one key feature: the values of the light sea-quark masses. Ensemble 1 is isospin symmetric, with the two light sea-quark masses equal and fixed to /2. Ensemble 2 features isospin breaking; here the two light-quark masses have the same average light-quark mass as ensemble 1, but the ratio of the light sea-quark masses is fixed to the value of determined from the MILC Collaboration’s study of pion and kaon electromagnetic mass splittings within the quenched approximation of QED Basak et al. (2016). Because the up and down sea-quark masses on this ensemble each take their physical values, a chiral extrapolation is not needed in our analysis. Comparing results on the two ensembles enables us to quantify the (tiny) effects of sea isospin breaking.

Our analysis strategy closely follows that of Ref. Chakraborty et al. (2017). On each ensemble, we calculate vector-current correlators with zero spatial momentum and all four combinations of local and spatially smeared interpolating operators at the source and sink. The smearing function is given in Eq. (A1) of Ref. Chakraborty et al. (2017), and we employ the same smearing parameters as in that work. To determine the quark-mass dependence of , we compute correlators with three valence-quark masses . With staggered quarks, the local vector current is not the conserved current, and must be renormalized. The renormalization factor for HISQ quarks, however, has only mild quark-mass dependence so it cancels when the strong-isospin-breaking correction is calculated as a percentage shift. Additional details on the correlator construction and wavefunction smearings can be found in Ref. Chakraborty et al. (2017).

We fit the matrix of correlators together using the multiexponential parametrization in Eq. (A6) of Ref. Chakraborty et al. (2017). The fit function includes contributions from both vector and opposite-parity states that arise with staggered valence quarks. The smeared correlators have smaller overlap with excited states than the local-local correlator, and therefore improve the determination of the energies and amplitudes. We fit the correlators over the symmetric time range [], thereby ensuring that the fit describes the correlator over the entire lattice time extent . To reduce the degrees of freedom in the fit, in practice we average the correlator at times and and fit only to the lattice midpoint; we also average the smeared-source, local-sink correlator with the local-source, smeared-sink correlator. Because our limited number of configurations do not enable us to reliably determine the smallest eigenvalues of the correlation matrix, we employ singular-value-decomposition (SVD) cuts with the values chosen to obtain stable fits with good correlated values. In practice, we replace all eigenvalues below the cut with the value of the SVD cut times the largest eigenvalue; this prescription increases the variance of the eigenmodes associated with the replaced eigenvalues and, thus, the errors on the fit parameters. We choose the number of states and fit range based on the stability of the ground-state and first-excited-state energies and amplitudes.

For both ensembles and all valence-quark masses, we obtain good correlated fits with stable central values and errors using , , and an SVD cut of 0.015, which modifies about 40% of the eigenvalues of the correlation matrix. For each of our six fits, the contribution to the from the 66 correlator data points ranges from about 45-80. Although the lowest-energy states in the vector-current correlators are pairs, we do not see any evidence of such states in our two-point correlator fits. This is not surprising because there are only a few states below the mass in these correlators, and their amplitudes are suppressed by the reciprocal of the spatial volume. The ground-state energies for the correlators with are MeV and MeV on the and ensembles, respectively; these are statistically consistent with the PDG average for the Breit-Wigner mass MeV Patrignani et al. (2016).

Following Ref. Chakraborty et al. (2017), we reduce the statistical errors in by replacing the correlator data at large times by the result of the multiexponential fit. Although the fit function is appropriate for the periodic lattice temporal boundary conditions, we correct for the finite lattice temporal size by using the infinite-time fit function and doubling the correlator extent to . We use the fitted correlator above fm; with this choice, roughly 80% of the value of comes from the data region. The values of computed with for fm agree within with those computed entirely from data, but with more than ten times smaller statistical errors for .

## Iii Analysis

We calculate using the method introduced by the HPQCD Collaboration Chakraborty et al. (2014), in which one constructs the and Padé approximants for the renormalized hadronic vacuum polarization function [] from time moments of zero-momentum vector-current correlation functions. These moments are proportional to the coefficients in a Taylor expansion of around . The true result is guaranteed to lie between the and Padé approximants. We employ the Padé approximant for obtained from the first six Taylor coefficients; the values of computed from the and Padé approximants differ by .

In Ref. Chakraborty et al. (2017), the and Padé approximants for are constructed from rescaled Taylor coefficients , where is the ground-state energy obtained from the two-point correlator fits. The rescaling was found to reduce the valence-quark-mass dependence of because the -meson pole dominates the vacuum polarization. In addition, the rescaling cancels most of the error from the uncertainty on the lattice scale , which enters via the muon mass present in the one-loop QED integral for . Figure 2 shows on -flavor ensemble at the up, down, and average light-quark masses. The valence-quark-mass dependence is statistically well resolved because the three points are strongly correlated, and is smaller after rescaling.

The physical value of the light-quark-connected contribution to is given by the sum of with two up quarks in the vector current and with two down quarks in the vector current weighted by the square of the quarks’ electromagnetic charges:

(1) |

We then define the absolute shift with respect to the isospin-symmetric value as , and the relative correction to be

(2) |

where .

Table 2 summarizes the isospin-breaking shifts on the and ensembles, both before and after rescaling the Taylor coefficients. As expected, we do not observe any significant difference between the two ensembles. The leading sea isospin-breaking contributions to are quadratic in the difference ; taking MeV gives a rough power-counting estimate of their size as . The differences in are smaller with rescaling because the valence-quark-mass dependence is milder.

The errors on the shifts in Table 2 stem primarily from the two-point correlator fits. The parametric errors from the lattice spacing are about a percent before rescaling, and are twenty times smaller with rescaling. The parametric errors from the current renormalization factor are 0.2%. The uncertainty due to the use of Padé approximants, which we take to be the difference between obtained with the and approximants, is about a percent. The 2.7% uncertainty on the ratio in Ref. Basak et al. (2016) stems largely from the estimate of electromagnetic effects, and leads to errors of about 2% and 1% on the physical up- and down-quark masses, respectively. Propagating the tuned quark-mass uncertainties to the physical using the measured slope of with respect to valence-quark mass changes the shifts in Table 2 by – () without (with) rescaling. Finally, the leading finite-volume and discretization effects, which arise from one-loop diagrams with intermediate states, cancel in because the charged pions in the loop are sensitive to the average of the up- and down-quark masses. Higher-order contributions are suppressed by , where is a typical chiral perturbation theory scale. We therefore estimate the systematic uncertainties in the shifts in Table 2 due to finite-volume and discretization effects to be 1% times the leading contributions, or .

Because the sea-quark-mass dependence of is tiny, we can compare the shift in the “direct” points in Fig. 2 to the valence-quark-mass dependence observed in Ref. Chakraborty et al. (2017), which analyzes several isospin-symmetric MILC HISQ ensembles at three lattice spacings and with a wide range of pion masses. Figure 3 of that work shows that the “raw” data for are approximately linear in from MeV down to the physical value with a slope that is independent of the lattice spacing. We can therefore estimate the change in that would result from varying between and from the unphysically heavy data in Ref. Chakraborty et al. (2017), and find a value consistent with the difference obtained from our fully physical calculation here.

The shifts in Table 2 only include contributions from quark-connected diagrams, with quark-disconnected contributions expected to be suppressed by . We estimate the quark-disconnected contribution to the strong-isospin-breaking correction from one-loop diagrams, which are especially sensitive to changes in the quark masses, within finite-volume chiral perturbation theory. Including the effect of taste splittings between the sea pions, which reduce the isopsin-breaking shift, we obtain for the -loop contribution . To account for resonance and higher-order contributions, we take about three times this value, or , as the uncertainty on the isospin-breaking shifts in Table 2 from missing quark-disconnected contributions. This conservative error estimate is approximately the size of the full quark-disconnected contribution to obtained by the BMW Collaboration on their coarsest ensemble with fm and similar taste splittings Borsanyi et al. (2017a); we expect the quark-disconnected contribution to the strong-isospin splitting to be smaller.

We obtain our final results for the relative correction by averaging the values on the two ensembles.

direct | with rescaling | |
---|---|---|

2+1+1 | +7.7(3.7) | +1.9(4.0) |

1+1+1+1 | +9.0(2.3) | +2.3(2.5) |

## Iv Result and outlook

We obtain for the relative strong isospin-breaking correction to the light-quark connected contribution to the muon hadronic vacuum polarization

direct, | (3) | ||||

with rescaling, | (4) |

where the errors include Monte Carlo statistics and all systematics.
Our result without rescaling the Taylor coefficients is consistent with phenomenological estimates of the dominant isospin-breaking contribution from – mixing using data Wolfe and Maltman (2011); Jegerlehner and Szafron (2011); Jegerlehner (2017b), –5, and chiral perturbation theory Cirigliano et al. (2002), , although – mixing will also include effects from quark-line disconnected diagrams that we
do not consider here.
Recent exploratory lattice-QCD calculations obtain somewhat smaller estimates for the relative strong isospin-breaking correction of roughly 0.2%–0.6% for MeV Giusti et al. (2017); Boyle et al. (2017).
We cannot directly compare our result in Eq. (3) with these values, however, because they were obtained with unphysically heavy pions and do not yet include systematic uncertainties.^{1}^{1}1Since our paper appeared, the RBC/UKQCD Collaboration obtained a new result for the strong-isospin-breaking shift at the physical pion mass of Blum et al. (2018), which agrees with our “direct” values in Table 2.

The percentage shifts in Eqs. (3) and (4) can be used to correct any existing or future result for the connected contribution to the hadronic vacuum polarization obtained with degenerate light quarks. Results for obtained without rescaling the Taylor coefficients should be corrected using Eq. (3); this applies to most recent lattice-QCD calculations. Equation (4) should be used to correct when rescaling is employed.

We have performed the first direct calculation of the strong-isosopin-breaking correction to at the physical up- and down-quark masses. We obtain an uncertainty on the relative correction of 0.7, which is smaller, and also more reliable, than the phenomenological estimate used in recent lattice-QCD calculations with equal up- and down-quark masses Chakraborty et al. (2017); Borsanyi et al. (2017b); Lellouch (2017). Thus, it reduces a significant source of uncertainty in , and is a crucial milestone towards a complete ab-initio lattice-QCD calculation of the hadronic contributions to with the sub-percent precision needed by the Muon and planned J-PARC experiments.

To improve our results in Eqs. (3) and (4), we will include quark-disconnected contributions, which are the dominant source of uncertainty, in a future work. We will also calculate directly the electromagnetic correction to using dynamical QCD+QED gauge-field configurations to be generated soon by the MILC Collaboration with quarks, gluons, and photons in the sea Zhou and Gottlieb (2014).

###### Acknowledgements.

We thank John Campbell, Vera Gülpers, Fred Jegerlehner, Laurent Lellouch, and Silvano Simula for useful discussions. Computations for this work were carried out with resources provided by the USQCD Collaboration, the National Energy Research Scientific Computing Center and the Argonne Leadership Computing Facility, which are funded by the Office of Science of the U.S. Department of Energy; and with resources provided by the National Institute for Computational Science and the Texas Advanced Computing Center, which are funded through the National Science Foundation’s Teragrid/XSEDE Program. Computations were also carried out on the Darwin Supercomputer at the DiRAC facility, which is jointly funded by the U.K. Science and Technology Facility Council, the U.K. Department for Business, Innovation and Skills, and the Universities of Cambridge and Glasgow. This work utilized the RMACC Summit supercomputer, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University. The Summit supercomputer is a joint effort of the University of Colorado Boulder and Colorado State University. This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. This work was supported in part by the U.S. Department of Energy under grants No. DE-AC05-06OR23177 (B.C.), No. DE-SC0010120 (S.G.), No. DE-SC0015655 (A.X.K.), No. DE-SC0009998 (J.L.), No. DE-SC0010005 (E.T.N.), No. DE-FG02-13ER41976 (D.T.), by the U.S. National Science Foundation under grants PHY14-17805 (J.L.), PHY14-14614 (C.D., A.V.), PHY13-16222 (G.P.L.), PHY12-12389 (Y.L.), and PHY13-16748 and PHY16-20625 (R.S.); by the Royal Society, STFC and Wolfson Foundation (C.T.H.D., D.H., J.K.); by the MINECO (Spain) under grants FPA2013-47836-C-1-P and FPA2016-78220-C3-3-P (E.G.); by the Junta de Andalucía (Spain) under grant No. FQM-101 (E.G.) by the Fermilab Distinguished Scholars Program (A.X.K.); by the German Excellence Initiative and the European Union Seventh Framework Program under grant agreement No. 291763 as well as the European Union’s Marie Curie COFUND program (A.S.K.); by the Blue Waters PAID program (Y.L.); and by the U.K. STFC under grants ST/N005872/1 and ST/P00055X/1 (C.M.). Brookhaven National Laboratory is supported by the U.S. Department of Energy under contract DE-SC0012704. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the United States Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.## References

- Bennett et al. (2006) G. W. Bennett et al. (Muon g-2), Phys. Rev. D73, 072003 (2006), arXiv:hep-ex/0602035 [hep-ex] .
- Blum et al. (2013) T. Blum, A. Denig, I. Logashenko, E. de Rafael, B. Lee Roberts, T. Teubner, and G. Venanzoni, (2013), arXiv:1311.2198 [hep-ph] .
- Grange et al. (2015) J. Grange et al. (Muon g-2), (2015), arXiv:1501.06858 [physics.ins-det] .
- Jegerlehner (2017a) F. Jegerlehner, (2017a), arXiv:1705.00263 [hep-ph] .
- Keshavarzi et al. (2017) A. Keshavarzi, D. Nomura, and T. Teubner, (2017), The KNT17 approach to calculating the HVP contribution to .
- Davier et al. (2017) M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, (2017), arXiv:1706.09436 [hep-ph] .
- Aubin and Blum (2007) C. Aubin and T. Blum, Phys. Rev. D75, 114502 (2007), arXiv:hep-lat/0608011 [hep-lat] .
- Chakraborty et al. (2017) B. Chakraborty, C. T. H. Davies, P. G. de Oliviera, J. Koponen, G. P. Lepage, and R. S. Van de Water, Phys. Rev. D96, 034516 (2017), arXiv:1601.03071 [hep-lat] .
- Della Morte et al. (2017) M. Della Morte, A. Francis, V. Gülpers, G. Herdoíza, G. von Hippel, H. Horch, B. Jäger, H. B. Meyer, A. Nyffeler, and H. Wittig, JHEP 10, 020 (2017), arXiv:1705.01775 [hep-lat] .
- Lellouch (2017) L. Lellouch (BMW), (2017), LO-HVP contribution to the muon from the Budapest-Marseille-Wuppertal collaboration.
- Lehner (2017) C. Lehner (2017) arXiv:1710.06874 [hep-lat] .
- Burger et al. (2014) F. Burger, X. Feng, G. Hotzel, K. Jansen, M. Petschlies, and D. B. Renner (ETM), JHEP 02, 099 (2014), arXiv:1308.4327 [hep-lat] .
- Wolfe and Maltman (2011) C. E. Wolfe and K. Maltman, Phys. Rev. D83, 077301 (2011), arXiv:1011.4511 [hep-ph] .
- Jegerlehner and Szafron (2011) F. Jegerlehner and R. Szafron, Eur. Phys. J. C71, 1632 (2011), arXiv:1101.2872 [hep-ph] .
- Jegerlehner (2017b) F. Jegerlehner, Springer Tracts Mod. Phys. 274, pp.1 (2017b).
- Hagiwara et al. (2004) K. Hagiwara, A. D. Martin, D. Nomura, and T. Teubner, Phys. Rev. D69, 093003 (2004), arXiv:hep-ph/0312250 [hep-ph] .
- Davier et al. (2010) M. Davier, A. Hoecker, G. Lopez Castro, B. Malaescu, X. H. Mo, G. Toledo Sanchez, P. Wang, C. Z. Yuan, and Z. Zhang, Eur. Phys. J. C66, 127 (2010), arXiv:0906.5443 [hep-ph] .
- Basak et al. (2016) S. Basak et al. (MILC), PoS LATTICE2015, 259 (2016), arXiv:1606.01228 [hep-lat] .
- Boyle et al. (2017) P. Boyle, V. Gülpers, J. Harrison, A. Jüttner, C. Lehner, A. Portelli, and C. T. Sachrajda, JHEP 09, 153 (2017), arXiv:1706.05293 [hep-lat] .
- Giusti et al. (2017) D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo, and S. Simula, in 35th International Symposium on Lattice Field Theory (Lattice 2017) Granada, Spain, June 18-24, 2017 (2017) arXiv:1710.06240 [hep-lat] .
- Bazavov et al. (2013) A. Bazavov et al. (MILC), Phys. Rev. D87, 054505 (2013), arXiv:1212.4768 [hep-lat] .
- Follana et al. (2007) E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), Phys. Rev. D75, 054502 (2007), arXiv:hep-lat/0610092 [hep-lat] .
- Borsanyi et al. (2012) S. Borsanyi et al., JHEP 09, 010 (2012), arXiv:1203.4469 [hep-lat] .
- Dowdall et al. (2013) R. J. Dowdall, C. T. H. Davies, G. P. Lepage, and C. McNeile (HPQCD), Phys. Rev. D88, 074504 (2013), arXiv:1303.1670 [hep-lat] .
- Patrignani et al. (2016) C. Patrignani et al. (Particle Data Group), Chin. Phys. C40, 100001 (2016).
- Chakraborty et al. (2014) B. Chakraborty, C. T. H. Davies, G. C. Donald, R. J. Dowdall, J. Koponen, G. P. Lepage, and T. Teubner (HPQCD), Phys. Rev. D89, 114501 (2014), arXiv:1403.1778 [hep-lat] .
- Borsanyi et al. (2017a) S. Borsanyi et al. (Budapest-Marseille-Wuppertal), (2017a), arXiv:1711.04980 [hep-lat] .
- Cirigliano et al. (2002) V. Cirigliano, G. Ecker, and H. Neufeld, JHEP 08, 002 (2002), arXiv:hep-ph/0207310 [hep-ph] .
- Blum et al. (2018) T. Blum, P. A. Boyle, V. Gülpers, T. Izubuchi, L. Jin, C. Jung, A. Jüttner, C. Lehner, A. Portelli, and J. T. Tsang, (2018), arXiv:1801.07224 [hep-lat] .
- Borsanyi et al. (2017b) S. Borsanyi, Z. Fodor, T. Kawanai, S. Krieg, L. Lellouch, R. Malak, K. Miura, K. K. Szabo, C. Torrero, and B. Toth, Phys. Rev. D96, 074507 (2017b), arXiv:1612.02364 [hep-lat] .
- Zhou and Gottlieb (2014) R. Zhou and S. Gottlieb (MILC), PoS LATTICE2014, 024 (2014), arXiv:1411.4115 [hep-lat] .